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Abstract. We present a study of correlations between D and D mesons produced in 
500 GeV/c 7r~-nucleon interactions, based on data from experiment E791 at Fermilab. 
We have fully reconstructed 791 ±44 charm meson pairs to study correlations between 
the transverse and longitudinal momenta of the two D mesons and the relative pro- 
duction rates for different types of D meson pairs. We see slight correlations between 
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the longitudinal momenta of the D and the D, and significant correlations between 
the azimuthal angle of the D and the D. The experimental distributions are compared 
to a next-to-leading-order QCD calculation and to predictions of the Pythia/Jetset 
Monte Carlo event generator. We observe less correlation between transverse momenta 
and different correlations between longitudinal momenta than these models predict 
for the default values of the model parameters. Better agreement between data and 
theory might be achieved by tuning the model parameters or by adding higher order 
perturbative terms, thus contributing to a better understanding of charm production. 

The relative production rates for the four sets of charm pairs, D°D°, D°D~ , D + D°, 
D + D~ , as calculated in the Pythia/Jetset event generator with the default parame- 
ters, agree with data as far as the relative ordering, but predict too many D + D° pairs 
and too few D + D~ pairs. 

PACS: 12.38.-t, 13.85.-t, 13.87.Ce, 13.87.Fh 



1 Introduction 

Using data from experiment E791 at Fermilab, we reconstruct pairs of charm 
mesons produced in 500 GeV/c 7r~-nucleon interactions, where y/s = 30.6 GeV, 
and use correlations between the mesons to probe two aspects of the hadropro- 
duction of mesons containing a heavy quark: the dynamics of the production 
of heavy quark-antiquark pairs and the subsequent hadronization of the quarks 
into hadrons. Correlations between the D and D momenta transverse to the 
beam direction are sensitive to corrections to the leading-order calculations of 
the cc cross section. Correlations between the longitudinal momenta, as well as 
differences in the production rates of the four types of DD pairs (D°D°, D°D~ , 
D + D°, and D + D~), provide information regarding the role of the remnants of 
the colliding hadrons in the hadronization process that transforms the charm 
quarks into charm mesons. 

In most studies of the hadroproduction of charm particles, distributions for 
single charm particles are used to probe the underlying production physics [1, 2]. 
The variables used to describe the single particle distributions are the transverse 
momentum with respect to the beam direction, p t , and either the rapidity y or 
the Feynman scaling variable Xf, where 

ii= — In [ — — Pz | and (1) 

x F = Pz /pT ax ^2 Pz /^- s . (2) 

E and p z are the center-of-mass energy and longitudinal momentum of the charm 
particle and y/s is the total center-of-mass energy. The center of mass is that of 
the pion-nucleon system. Such single charm studies are insensitive to correlations 
between the two charm hadrons in a single event. 

We have fully reconstructed 791 ± 44 DD pairs. Based on this sample, we 
present background-subtracted, acceptance-corrected distributions for the fol- 
lowing variables: 

1. the invariant mass of the pair of charm mesons, M D -^; 
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2. the square of the vector sum of the transverse momenta, with respect to 
the beam direction, of the D and D mesons (p^ D — = \pt,D + P t jj\ 2 )', 

3. correlations between xf,d and x F p, as well as yu and y-^; 

4. Axp = x pd — x F p and Sxp = xp.D + x fd'i 

5. Ay = y D - y^ and s V = Vd + J/jj; 

6. correlations between the squares of the magnitudes of the transverse mo- 
menta of the D and D mesons, p\ D and P^^] 

7- Ap\ = \pl D - p\-\ and Ep* = p^+p*-; 

8. the azimuthal separation between the momentum vectors of the D and D 
mesons in the plane perpendicular to the beam direction, A(f> = (minimum 
of \(f> D - <%| and 360° - \<j> D - <%|); 

9. correlations between the azimuthal separation (A(j)) and the scalar sum 
and difference of the D and D transverse momenta, Ap\ and Ylp\ ; 

In addition, this paper reports the relative production rates for each type 
of DD pair (D°D a , D°D~, D+D°, and D+D~), and compares the rapidity 
correlations for the various DD pair combinations. 

We also investigate the extent to which the observed charm-pair correlations 
can be duplicated by simply convoluting the observed single charm particle dis- 
tributions. In addition, we compare our measured distributions to three sets of 
theoretical predictions: 

1. the distributions of cc pairs from a next-to- leading-order perturbative QCD 
calculation by Mangano, Nason and Ridolfi[3, 4]; 

2. the distributions of cc pairs from the Pythia/Jetset Monte Carlo event 
generator [5] which uses a parton-shower model to include higher-order per- 
turbative effects [6]; and 

3. the distributions of DD pairs from Pythia/Jetset which uses the Lund 
string model to transform cc pairs to DD pairs [7]. 

In Table 1, we compare the E791 charm-pair sample to those from other fixed- 
target experiments (both hadroproduction and photoproduction). The largest 
previous sample of fully-reconstructed hadroproduced charm pairs used to study 
correlations is 20 pairs from the CERN 7r~-nucleon experiment NA32[10]. Some 
studies have been conducted with partially-reconstructed charm hadrons, in 
which the direction but not necessarily the magnitude of the charm particle 
momentum is determined directly. NA32 partially reconstructed 642 such charm 
pairs [11]. In photoproduction experiments, the largest sample of charm pairs 
reconstructed is from the E687 data[15], with 325 fully-reconstructed and 4534 
partially-reconstructed charm pairs. In the E687 partially-reconstructed sample, 
one D meson is fully reconstructed and the momentum vector of the other charm 
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Table 1. Summary of fully-reconstructed and partially- reconstructed charm-pair 
samples from hadroproduction and photoproduction fixed-target experiments. In 
the last column, we list the physics variables studied in each experiment. The 
variables are defined in the text. 



Experiment 


Beam 


Number 


Measured Pair 




Energy (GeV), 


of Pairs 


Variables 




Beam Type, 


Reconstructed 






and Target 






E791 


500 7T" 


791 fully 


A<t>, \Ep t \\ Spl 








Ap\, correlations, 


(this paper) 


Pt, C 




Exf, Axf, Ey, 








A V> M DW 








a D°D° : a D'>D- ■ 








a n+n° : a D+D- 


WA92 [8] 


350 7T- 


475 partially 1 


AcP, \Eptf, M Dl5 , 




Cu 




Axf, Exf, Ay 


E653 [9] 


800 p 


35 partially 


At, iA>tr, 




emulsion 




Ay, M D p, cos 6 cm 


NA32 [10] [11] 


230 7T- 


20 fully 


\Ep t \ 2 , Ay, M D -jj 


(ACCMOR) 


Cu 


642 partially 


\Ep t f, Ay, M n75 , A<f> 


WA75 [12] 


350 7T- 


177 partially 


A<t>, Ay 




cmulsion 


120 partially 


M D jj, Ex F , Epf 


NA27 [13] 


400 p 


17 fully 


\Ep t f, Ex F , Ay, M Dl5 


(LEBC) 


H 2 


107 partially 


A(j>, a D o D a : 








<T D»D-+D + D° 1 a D+D- 


NA27 [14] 


360 7T- 


12 fully 


\E Pt \ 2 , Ex F , 








Ay, M Dl5 


(LEBC) 


H 2 


53 partially 


A<p 


E687 [15] 


200 7 


325 fully 


Acf>, \Ep t \-', Ay, M Dl5 




Be 


4534 partially 


A<p, Ay, M D jj 


NA14/2 [16] 


100 7 


22 fully 


A<j>, \Ep t \, Ep z , 




Si 




Ay, m dd 
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meson is determined by scaling the momentum vector of low-momentum charged 
pions from the decays D* ± — ► D°n ± . 

In the analysis presented here, we have completed an extensive study of accep- 
tance corrections. Acceptance corrections are made as a function of the eight vari- 
ables that describe the D and D degrees of freedom: ((y,pt, <j>,n)£>, (y,p t ,(fi, n)-jj). 
Here n is the number of decay tracks from the D meson. Corrections are also 
made for the branching fractions of the reconstructed D and D decay modes. 

We performed a maximum likelihood fit to the two-dimensional reconstructed 
candidate D mass distribution, including terms in the likelihood function for the 
true DD pairs that are the signal of interest, and also terms for combinations 
of a true D with background, combinations of a true D with background, and 
combinations of two background candidates in the same event. From the full data 
set, the resulting number of true fully reconstructed DD pairs was 791 ± 44. In 
making the distributions for single-charm and charm-pair physics variables, a 
likelihood fit was performed for each bin of the relevant physics variable. 

Since we fully reconstruct both the D and D meson, our results have fewer 
systematic errors than previously published results based on partially recon- 
structed pairs. In particular, we do not need to correct for missing tracks or 
possible contamination from baryons. 

In the next section, we review the current theoretical understanding of the 
hadroproduction and hadronization of charm quarks. In the Appendix, we use 
both theoretical calculations and phcnomcnological models to investigate the de- 
pendence of various measurable properties of charm production on higher-order 
QCD effects, the charm quark mass, the parton distribution functions, the fac- 
torization scale and the renormalization scale. In Sec. 3, we describe the E791 
detector and data processing. In Sec. 4, we describe the optimization of selection 
criteria for charm pairs. We discuss the extraction of background-subtracted dis- 
tributions and corrections for acceptance effects in Sec. 5. In Sec. 6, we present 
the measured distributions for the charm pairs and compare them to the distri- 
butions predicted by (uncorrelated) single-charm distributions and to theoretical 
predictions. We summarize our results in Sec. 7. 

2 Theoretical Overview 

The charm quark is the lightest of the heavy quarks. Its relatively small mass 
ensures copious charm particle production at energies typical of fixed-target 
hadroproduction experiments. Its relatively large mass allows calculation of the 
large-momentum-transfer processes responsible for producing cc pairs using per- 
turbative quantum chromodynamics (QCD). The consequence of the charm 
quark being the lightest heavy quark — more specifically, having a mass not 
sufficiently larger than Aqcd — is that there are considerable uncertainties as- 
sociated with these calculations. Such large theoretical uncertainties, combined 
with conflicting experimental results from early charm hadroproduction experi- 
ments, have made systematic comparisons between theory and data difficult to 
interpret. Recent calculations of the full next-to-lcading-order (NLO) differential 
cross sections by Mangano, Nason and Ridolfi (MNR)[3] and others, as well as 
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unprecedented numbers of charm particles reconstructed by current fixed-target 
experiments, have allowed more progress to be made in this field. 

In this section, we outline the theoretical framework used to describe the 
hadroproduction of charm pairs, focusing on the framework used by the fol- 
lowing two packages: the FORTRAN program HVQMNR[4], which implements 
the MNR NLO perturbative QCD calculation for charm quarks, and the Py- 
thia/ Jetset Monte Carlo event generator [5], which makes predictions for charm 
particles based on leading order parton matrix elements, parton showers and the 
Lund string fragmentation model. In the Appendix we examine predictions from 
these two packages for the same beam type and energy as E79I for a wide range 
of theoretical assumptions to determine how sensitive the theoretical predictions 
arc to 

1. the inclusion of higher order terms (a^ or parton shower contributions); 
and 

2. non-perturbative effects, including 

(a) variations in parameters such as the mass of the charm quark and 
alternative parton distribution functions; 

(b) changes in the factorization and renormalization scales; and 

(c) other non-perturbative effects (hadronization and intrinsic transverse 
momentum of the colliding partons). 

2.1 Charm Quark Production 

Both the HVQMNR and Pythia/ Jetset packages use a perturbative QCD 
framework to obtain the differential cross section for producing a cc pair: 

dcr c - = ^2 dx b dx t tf(x b , li f ) fifa, li f ) da i:j (x b P b ,x t P t ,p c ,Pc,m c , n R ), (3) 
where 

• P b (P t ) is the momentum of the beam (target) hadron in the center of 
mass of the colliding hadrons; 

• x b (xt) is the fraction of P b (P t ) carried by the hard-scattering parton from 
the beam (target) hadron; 

• fi {ft) i s the parton distribution function for the beam (target) hadron; 

• fj,R, the renormalization scale, and /zp, the factorization scale, come from 
the perturbative QCD renormalization procedure which transforms the 
QCD coupling constant g — \/4ira s and the 7r~ and nucleon wave functions 
from "bare" (infinite) values to physical (i.e., finite and measurable) values; 

• d&ij is the differential cross section for two hard-scattering partons to pro- 
duce a pair of charm quarks, each with mass m c , and with four- momenta 
p c and pc- 

http : //link. springer .de/link/ service/ journals/10105/ index.html 



EPJdiroct CI (2008) Springer- Verlag 



73 



Leading order (a 2 s ) contributions to the cc cross section require the charm 
and anticharm quarks to be produced back-to-back in the center of mass of 
the cc pair. The (unknown) partonic center of mass is boosted in the beam di- 
rection with respect to the (known) hadronic center of mass. At fixed target 
energies, this boost smears the longitudinal momentum correlation while pre- 
serving the transverse correlation. Therefore, leading order calculations predict 
delta function distributions (i.e., maximal correlations) for variables which mea- 
sure transverse correlations, such as A<p c - — 180° and p% c - = 0, but predict small 
correlations in the longitudinal- momentum correlation variables Axf, Sxf, Ay 
and Ey. 

These leading-order predictions are altered by the inclusion of higher order 
effects. The HVQMNR program adds the NLO (af) corrections to the leading 
order calculation. NLO processes such as gg — ► ccg produce cc pairs that are no 
longer back-to-back, smearing the leading order delta function distributions for 
Acj) cZ and p t 2 c -. 

The Pythia/ Jetset event generator accounts for higher order perturba- 
tive QCD effects via a "parton shower" approach [17]. In this approach each of 
the two incoming and two outgoing partons, whose distributions are based on 
leading-order matrix elements, can branch — backwards and forwards in time 
respectively — into two partons, each of which can branch into two more par- 
tons, etc. This evolution continues until a small momentum scale is reached. In 
addition, the Pythia/ Jetset event generator gives the hard-scattering partons 
an intrinsic transverse momentum k t ■ Both the parton showers and the intrinsic 
transverse momentum tend to smear the transverse correlations, as shown in the 
Appendix. 

The extent to which transverse-momentum correlations are smeared provides 
a measure of the importance of higher order perturbative effects. In addition, 
since the leading order calculation predicts very little longitudinal-momentum 
correlation, an enhancement of the longitudinal-momentum correlation also pro- 
vides evidence for higher order perturbative effects or non-perturbative effects 
such as hadronization, described below. 

2.2 Hadronization 

The process whereby charm quarks are converted to hadrons is known as hadroniza- 
tion or fragmentation. Since this process occurs at an energy scale too low to 
be calculable by perturbative QCD, fragmentation functions are used to param- 
eterize the hadronization of the charm quark. Such functions have been mea- 
sured by several e + e _ experiments. The hadroproduction environment in 7r _ -N 
interactions, however, is quite different from the e + e~ environment. In e + e~ in- 
teractions, the light quarks in the produced charm hadrons must come from the 
vacuum. Hadroproduction leaves light quark beam and target remnants which 
are tied by the strong force to the charm quarks. Interactions between these rem- 
nants and the charm quarks can dramatically affect the momentum and flavor 
of the observed charm hadrons. 

The Pythia/ Jetset event generator uses the Lund string fragmentation 
framework, described in the Pythia/ Jetset manual [5], to hadronize the charm 
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quarks. To illustrate this model we consider an example from E791 where a gluon 
from a tt~ and a gluon from a nuclcon in the target interact to form a cc pair. 
This accounts for ~ 90% of the theoretical cross section for 500 GeV/c 7r~-N 
interactions. After the gluon-gluon fusion, the remnant 7r~ and nucleon are no 
longer color-singlet particles. The remnant n~ is split into two valence quarks, 
and the remnant nucleon into a valence quark plus a diquark. Given this minimal 
set of partons — (c, c), (u, d) n , and (qq, g)jv — the two dominant ways to make 
color-singlet strings, and the ones PYTHIA uses are [18]: 

(c, u n ), (c, gjv), and (d n , qqN); or (4) 
(c, d*), (c, qq N ), and (u n , q N )- 

In the center of mass of a particular qq system, such as cd, the c and d are 
moving apart along the string axis. As they move apart, energy is transferred 
to the color field. When this energy is great enough, qq pairs are created from 
the vacuum with equal and opposite transverse momentum (with respect to the 
string axis) according to a Gaussian distribution. The transverse momentum 
relative to the string axis of the resulting cq meson is determined by the q quark 
since the c contributes none. The longitudinal momentum of the meson is given 
by a fragmentation function which describes the probability that a meson will 
carry off a fraction z of the available longitudinal momentum. By default, heavy 
quark fragmentation is performed according to a Lund fragmentation function [7] 
modified by Bowler [19]: 

f(z) oc ^^expj— ±) (5) 

where m 2 = M% + p 2 is the transverse mass of the hadron and vhq is the 
mass of the heavy quark. The default Pythia/ Jetset settings are a = 0.3 and 
6 = 0.58(GeV/c 2 )- 2 . 

When the remaining energy in the string drops below a certain cutoff (depen- 
dent on the mass of the remaining quarks) a coalescence procedure is followed, 
which collapses the last partons into a hadron while conserving energy. The en- 
tire string system is then boosted back into the lab frame. In the case of a (c, d v ) 
or (c, 77^) string, this boost will tend to increase the longitudinal momentum of 
the charm hadron with respect to the charm quark since the d v and u n will tend 
to have large longitudinal momentum. The opposite will occur, however, for a 
(c, q N ) string. 

In some fraction of events, strings will be formed with too little energy to 
generate qq pairs from the vacuum. In these cases the c quark (antiquark) will 
coalesce into a single meson with the beam antiquark (quark) or will coalesce 
into a single baryon (meson) with the target diquark ("bachelor" quark). This 
will tend to enhance production of charm hadrons with a light quark in common 
with a valence beam quark in the forward direction (beam fragmentation region) 
and production of charm hadrons with a light valence quark or diquark in com- 
mon with the target in the target fragmentation region. This phenomenon has 
been used to explain the leading particle effect seen in charm hadroproduction 
experiments [20]-[25]. 
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Fig. 1. Scatter plots of jfc vs. y c and vs - Vd, from 100,000 Pythia/ Jetset 
cc and DI> events, showing the correlation introduced by the hadronization 
model. As discussed in Sec. 4, we only reconstruct DD events in the region 
-0.5 < y D ^ < 2.5. 



However, in most events, the string has sufficient energy 2 to produce at 
least one qq pair from the vacuum. In this type of beam/target "dragging," 
the strength of the dragging is not dependent on the light quark content of the 
produced particle. 

These effects are evident in Fig. 1, which shows a scatter plot of the charm 
and anticharm rapidities for Pythia/ Jetset DD events. 3 Comparison of the 
scatter plot of the charm and anticharm quark rapidities, Fig. la, to the scatter 
plot of the D and D rapidities, Fig. lb, clearly demonstrates that significant 
correlations are introduced by hadronization. 

Both the degree of correlation between the D and D longitudinal momenta as 
well as differences in production of the four types of DD pairs — D° D° , D° D~ , 
D + D° and D + D~ — provide information about the charm quark hadronization 
process in a hadronic environment. 



3 Experiment E791 

The results reported in this paper are based on a data sample recorded by Fermi- 
lab experiment E791 during the 1991/92 fixed-target run. The E791 spectrometer 
is illustrated in Fig. 2. A 500 GeV/c ir~ beam impinged on platinum and carbon 
targets. The spectrometer consisted of proportional wire chambers (PWC's) and 
silicon microstrip detectors (SMD's) upstream and downstream of the targets, 
two magnets, 35 drift chamber (DC) planes, two Cerenkov counters, an electro- 
magnetic calorimeter, a hadronic calorimeter and a muon detector composed of 
an iron shield and two planes of scintillation counters. 

2 In contrast, at high xp most of the particle energy is taken up by the individual partons, 
so that the string has insufficient energy to produce qq pairs, and large asymmetries are seen 
by experiments. 

3 Dcfault values are used for all Pythia/ Jetset parameters. 
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TAGGED PHOTON SPECTROMETER 
E791 




Fig. 2. The E791 spectrometer. 



The spectrometer was an upgraded version of the apparatus used in Fermilab 
experiments E516, E691, and E769[26]. The major differences between the earlier 
versions and E791 were the addition of more planes of SMD's, enhancement of 
the muon identification system, new front-end detector-signal digitizers and a 
new data acquisition system. In general, E791 finds that the most important 
parts of the spectrometer for analysis are the charged-particle tracking system 
and the threshold Cerenkov counters; although the threshold Cerenkov counters 
were used minimally in the analysis reported in this paper. 

3.1 Target 

The target consisted of five foils with center-to-center separations that varied 
from 14.8 to 15.4 mm. The most upstream foil was 0.5 mm thick and was made 
of platinum to provide a significant interaction probability in a thin target. The 
next four foils were 1.6 mm thick and were made of industrial diamond. The low 
Z of these carbon targets minimized multiple scattering, while the higher den- 
sity of diamond permitted thinner downstream targets for the same interaction 
probability. The total pion interaction length of all five targets was about 1.9%. 
This target arrangement was chosen so that most of the particles with lifetimes 
and momenta within the range of interest to this experiment have a decay vertex 
in the gaps, where there is less background from secondary interactions. 

3.2 The Spectrometer 

The 7T~ beam particle was tracked with eight PWC planes and six SMD planes 
upstream of the target region. Downstream of the targets, the charged-particle 
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tracking system consisted of 17 SMD planes, two PWC planes, and 35 drift 
chamber planes. In general several planes of tracking chambers with different 
angular orientations around the beam axis were grouped together in each track- 
ing station to provide hit ambiguity resolution. The various coordinates (x, y, 
w, u, v) measured by the planes in the tracking chamber stations were defined 
relative to a right-handed coordinate system x — y — z in which increasing z was 
in the beam direction, x was the horizontal dimension and y increases vertically 
upward. The w, u and v axes were rotated by +60°, +20.5°, and —20.5° with 
respect to the positive x axis. 

The beam PWC's[27] had a wire spacing of 1 mm and were arranged in two 
stations widely separated in z to measure the angle of the incoming beam particle 
with high precision. The first station was 31 m upstream, and the second was 
12 m upstream of the last carbon target. Each station consisted of 4 planes: two 
staggered x planes, a y plane and a w plane. 

The beam SMD's had a pitch of 25 /im and were also arranged in two stations, 
each with an x, y and w plane. The first SMD station was 80 cm upstream 
of the most downstream target, and the second station was 30 cm upstream 
of this target. The system of SMD's downstream of the targets started 2.8 cm 
downstream of the last target and extended for 45 cm. It had a maximum angular 
acceptance of about ±125 mr in both x and y. Each of the first two planes (x 
and y) had an active area of 2.5 cm by 5 cm, a pitch of 25 /um in the central 
9.6 mm and 50 fxm in the outer regions, and an efficiency of about 84%. The 
next nine planes were identical to those used in E691[28]. Each plane had a pitch 
of 50 /urn, and an efficiency from 88% to 98%. They were instrumented to give 
an acceptance of ±100 mrad with respect to the center of the most downstream 
target. They measured x — y — v — y — x — v — x — y — v coordinates respectively. 
The final six SMD planes had active areas of 9 cm by 9 cm. The inner 3 cm had 
a pitch of 50 /xm while the outer regions had an effective pitch of 200 /im. These 
measured v — x — y — x — y — v coordinates. The efficiencies ranged from 96% to 
99%. 

The drift chambers were arranged in four stations as illustrated in Fig. 2. 
Each station was subdivided into substations with plane orientations such that 
an x~ y — z space point could be reconstructed in each substation. The charac- 
teristics of these chambers are given in Table 2. Since the beam, which operated 
at about 2 MHz throughout the run, passed through the center of the drift cham- 
bers in a small region instrumented with very few wires, each plane had a central 
inefficient region in which the efficiency decreased to < 10% and the resolution 
was degraded by as much as a factor of four. The profile of the efficiency and 
resolution degradation region was approximately gaussian in an angular region 
of three to four mrad centered on the beam. The extent of the inefficient region 
increased with time during the run and is the major source of systematic uncer- 
tainty associated with the acceptance at large xp. Each substation of the first 
drift chamber station was augmented by a PWC which measured the y coordi- 
nate. These PWC's had a wire spacing of 2 mm. Typical inclusive single charm 
acceptances for two, three, and four particle D decays are shown in Fig. 3. 

Momentum analysis was provided by two dipole magnets that bent particles 
in the same direction in the horizontal plane. The transverse momentum kicks 
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Table 2. Characteristics of the DC tracking chambers for E791. "View" refers 
to the coordinate measured by that plane, with u = +20.5°, v = —20.5° and x' 
staggered by one-half a wire spacing relative to the x plane. The efficiency is for 
the region outside the central inefficient area. 



Station 


Dl 


D2 


D3 


D4 


Approximate size (cm) 


130 x 75 


280 x 140 


320 x 140 


500 x 250 


Number of substations 


2 


4 


4 


1 


Views per substation 


x, x', u, v 


u, x, v 


U, X, V 


u, x, v 


u and v cell size (cm) 


0.446 


0.892 


1.487 


2.974 


x cell size (cm) 


0.476 


0.953 


1.588 


3.175 


z position of first plane 


142.4 


381.4 


928.1 


1738. 


z position of last plane 


183.7 


500.8 


1047.1 


1749.2 


Approx. resolution (/xm) 


430 


320 


260 


500 


Typical efficiency 


92% 


93% 


93% 


85% 




Fig. 3. E791 acceptance functions vs. xp and pi for Kir, Kirn and Ktttttt 
candidates. The inclusive charm acceptance shown here was used to obtain 
approximately 150,000 reconstructed D° and D + charm meson decays in the 
E791 spectrometer as described, for example, in the study of single-charm pro- 
duction [29] and [30]. The p\ acceptance is obtained for charm mesons with 
—0.1 < xp < 0.6. The acceptance for events in which both charm particles are 
detected is quite different and is documented in detail later in this paper. (See 
Section 5.2.) 
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were 212 MeV/c for the first magnet and 320 MeV/c for the second magnet. 
The centers of the two magnets were 2.8 m and 6.2 m downstream of the last 
target, respectively. The x — y — z aperture of the pole faces of the first magnet 
was 183 cm by 81 cm by 100 cm and that of the second magnet was 183 cm by 
86 cm by 100 cm. 

Two segmented, gas-filled, threshold Cerenkov counters[31] provided particle 
identification over a large range of momenta. The threshold momenta above 
which a charged particle emits light were 6, 20 and 38 GeV/c for 7r's, K's, and 
p's, respectively, for the first counter, and 11, 36, and 69 GeV/c for the second. 
The particle identification algorithm correlates the Cerenkov light observed in a 
given mirror-phototube segment with the charged particle tracking information. 
The algorithm indicates the likelihood that a charged particle of a given mass 
could have generated the observed Cerenkov light in the segment (s) in question. 

The electromagnetic calorimeter, which we called the Segmented Liquid Ion- 
ization Calorimeter (SLIC) , consisted of 20 radiation lengths of lead and liquid 
scintillator and was 19 m from the target. Layers of scintillator counters 3.17 and 
6.24 cm wide were arranged transverse to the beam and their orientations al- 
ternated among horizontal and ±20.5° with respect to the vertical direction[32]. 
The hadronic calorimeter consisted of six interaction lengths of steel and acrylic 
scintillator. There were 36 layers, each with a 2.5-cm-thick steel plate followed by 
a plane of 14.3-cm-wide by 1-cm-thick scintillator slats; the slats were arranged 
alternately in the horizontal and vertical directions, and the upstream and down- 
stream halves of the calorimeter were summed separately [33]. The signals from 
the hadronic calorimeter as well as those from the electromagnetic calorimeter 
were used for electron identification. Signals from both calorimeters were used 
to form the transverse energy requirement in the hardware trigger [34] . Electron 
identification was not used in this analysis. 

Muons were identified by two planes of scintillation counters located behind 
a total of 15 interaction lengths of shielding, including the calorimeters. The first 
plane, 22.4 m from the target, consisted of twelve 40-cm-wide by 300-cm-long 
vertical scintillation counters in the outer region and three counters 60 cm wide in 
the central region. The second plane, added for E791, consisted of 16 scintillation 
counters 24.2 m from the target. These counters were each 14 cm wide and 
300 cm long, and measured position in the vertical plane. These counters were 
equipped with TDC's which provided information on the horizontal position of 
the incident muons [34]. Muon identification was not used in this analysis. 

3.3 Trigger and Data Acquisition 

To minimize biasing the charm data sample, the trigger requirements were very 
loose. The most significant requirements were that the signal in a scintillation 
counter downstream of the target be at least four times the most likely signal 
from one minimum-ionizing particle, and that the sum of the energy deposited in 
the electromagnetic and hadronic calorimeters, weighted by the sine of the angle 
relative to the beam, be above a threshold corresponding to 3 GeV of transverse 
energy. The time for the full hardware trigger decision was about 470 ns. This 
trigger was fully 100% efficient for charm decays. 
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A total of 24,000 channels were digitized and read out in 50 /us with a parallel- 
architecture data acquisition system[35]. Events were accepted at a rate of 9 kHz 
during the 23-second Tevatron beam spill. The typical recorded event size was 
2.5 kbytes. Data were written continuously (during both the 23-second spill and 
the 34-second interspill periods) to forty- two Exabyte [36] model 8200 8-mm- 
tape drives at a rate of 9.6 Mbytes/s. Over 2 x 10 10 hadronic interactions were 
recorded on 24,000 tapes. 

3.4 Data Processing 

The 2 x 10 10 interactions recorded constitute about 50 Terabytes of data. Event 
reconstruction and filtering took place over a period of two and a half years at 
four locations: the University of Mississippi, The Ohio State University (moved 
to Kansas State University in 1993), Fermi National Accelerator Laboratory, 
and O Centro Brasileiro de Pesquisas Fi'sicas, Rio de Janeiro (CBPF). The first 
three sites used clusters of commercial UNIX/RISC workstations controlled from 
a single processor with multiprocessor management software [3 7], while CBPF 
used ACP-II custom-built single-board computers [38]. 

As part of the reconstruction stage, a filter was applied which kept ~20% 
of the events. This filter was effectively an offline trigger. To pass this filter, an 
event was required to have a reconstructed primary production vertex whose 
location coincided with one of the target foils. The event also had to include at 
least one of the following: 

1. At least one reconstructed secondary decay vertex of net charge for an 
even number of decay tracks and ±1 for an odd number of decay tracks. 
The longitudinal separation of the secondary vertex from the primary had 
to be at least four sigma for secondary vertices with three or more tracks 
and at least six sigma for vertices with two tracks, where sigma is the error 
in the separation, 

2. At least one reconstructed K s — ► n~n + or A — > pit candidate whose decay 
was observed upstream of the first magnet, 

3. and for part of the run, at least one reconstructed <j> — ► K + K~ candidate. 

For one-third of the data sample, several other classes of events were also kept. 
These are included in analyses not covered in this paper: 

4. Events in which the net charge of all the reconstructed tracks was negative 
and their total momentum was a large fraction of the beam momentum. 

5. K s — > 7r~7r + or A — > pir candidates that decayed inside the aperture of the 
first magnet. 

Following the initial reconstruction/filter, which was applied to all events, ad- 
ditional selections of events were made to further divide the large data sample 
into subsets by class of physics analysis. 
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3.5 Detector Performance 

The important detector performance characteristics for this analysis are the 
resolution for reconstructing the positions of both the primary interaction and 
secondary decay vertices, the efficiency for reconstructing the trajectories of 
charged particles, the resolution for measuring charged track momenta, and the 
efficiency and the misidentffication rates for identifying charged pions and kaons 
using information from the Cerenkov counters. 

The resolution for measuring the position of the primary vertex along the 
beam direction varies from about 240 fim for the most downstream target foil 
to 450 /im for the upstream foil. The variation is due to the extrapolation from 
the SMD system and to multiple scattering in material downstream of the inter- 
action. The mean number of reconstructed tracks used to fit the primary vertex 
is seven. The measured secondary vertex resolution depends on the decay mode, 
the momentum of the D, and the selection criteria. For example, the vertex 
resolutions along the beam direction for K~n + and K~Tr + n~7r + arc 320 and 
395 /im, respectively, for a mean D° momentum of 65 GeV/c, and worsen by 33 
and 36 fim for every 10 GeV/c increase in D° momentum. 

The total efficiency, including acceptance, for reconstructing charged tracks 
is approximately 80% for particles with a momentum greater than 30 GeV/c and 
drops to around 60% for particles of momentum 10 GeV/c. (This includes the 
inefficiency in the beam region.) For tracks which pass through both magnets and 
have a momentum greater than 10 GeV/c, the average resolution for measuring 
charged particle momentum p is Sp/p = 0.6% © (0.02p)% where stands for 
the quadratic sum, and p is in GeV/c. Tracks which pass through only the first 
magnet have a resolution Sp/p — 2% © (0.1p)%. The mean D mass resolution 
for hadronic decays to two, three and four charged particles varies from 13 to 
8 MeV/c 2 as the decay multiplicity increases. The mass resolution varies by 
about a factor of 2 between low and high momentum D mesons. 

In most E791 analyses, the Cerenkov counters play a very important role 
[39]. However, in this analysis with the two fully reconstructed D-meson decays, 
the Cerenkov counters play a minimal role. We use the Cerenkov counters for 
charged kaon identification. The kaon identification efficiencies and misidentffica- 
tion probabilities vary with longitudinal and transverse momentum and with the 
signatures required in the Cerenkov counters. For typical particle momenta in 
the range 20 GeV/c to 40 GeV/c, and for the nonstringent criteria used for some 
of the final states, in this analysis, the Cerenkov identification efficiency for a 
kaon ranged from 64% to 72% while the probability for a pion to be misidentificd 
as a kaon ranged from 6% to 12%. 

A complete Monte Carlo simulation of the apparatus was used in this analysis 
to calculate the efficiency and investigate systematic effects. The simulation in- 
cluded all relevant physical processes such as multiple interactions and multiple 
scattering as well as geometrical apertures and resolution effects. It produced 
data in the same format as the real experiment. That Monte Carlo data was 
then reconstructed and analyzed with the same software as the real data. 
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4 Event Selection 

In each E791 event, we search for two charm mesons (D°, D°, D + or D~) de- 
caying to Cabibbo-favored final states that can be reconstructed with relatively 
high efficiency: D° -> K~tt + , D° -> K~tz + 'k-'k + , D+ -> K-ir+ir+, and the 
charge conjugate modes. 4 To optimize the efficiency for reconstructing charm 
pairs, we search for both D candidates simultaneously, rather than searching for 
the two candidates consecutively. In such a simultaneous search, we can require 
that one candidate or the other satisfy a fairly stringent selection criterion based 
on a particular variable used to discriminate charm decays from background, or 
that both candidates satisfy less stringent criteria. 

We start with a sample of events, each containing two candidate K~tt + , 
K~it + it + or K~it + tt~it + combinations with invariant mass between 1.7 and 
2.0 GeV/c 2 , and rapidity in the range —0.5 < y D jj < 2.5. These candidates are 
found by looping over all reconstructed tracks. The primary vertex is then refit 
after removing tracks which are associated with either candidate. No particle 
identification requirements are applied at this time. Candidates are rejected if 
any charged track, the primary vertex or either of the two secondary vertices do 
not meet minimal fit quality criteria. The sample of candidate pairs that pass 
just these criteria is dominated by combinatoric backgrounds. To choose further 
selection criteria, we use this sample to represent background. To represent sig- 
nal, we use reconstructed charm pairs generated with the Monte Carlo program 
described at the end of the previous section. We then search for selection criteria 
that provide optimal discrimination between signal and background. 

In order to extract the signal, we use selection criteria defined by discrimi- 
nation variables (properties of the candidate event) and minimum or maximum 
allowed values for each variable. For candidate pairs with the same final states 
(i.e., both Ktt, both Kirir, or both Kiririr), the same discrimination variables 
and maximum or minimum values are used for both candidate D's; for pairs with 
different final states, the discrimination variables are allowed to be different for 
the two candidate D's. 

The discrimination variables used address the following questions. Is a D 
candidate consistent with originating from the primary interaction vertex? Is 
the vertex formed by the decay products of a D candidate well separated from 
the primary interaction vertex and not inside a target foil? Do any of the decay 
products of the D candidate appear to originate from the primary interaction 
vertex or from the other D candidate vertex? Is the scalar sum of the squares 
of the transverse momenta of the D candidate decay products, with respect to 
the D candidate trajectory, indicative of a heavy meson decay? Is the Cerenkov 
information for the kaon candidate consistent with that for a real kaon? As an 
example of the cuts used: the p t balance cut was 400 McV/c, and the secondary 
vertices were separated from the primary vertices by 8 times the rms- uncertainty 
in the separation. 

To optimize the significance of the signal, we repeatedly choose the selection 
criterion that maximizes Ns/^Ns + Nb while rejecting no more than 5% of the 
(Monte Carlo) signal. N$ is the number of signal pairs satisfying the selection 

4 Unless noted otherwise, charge conjugate modes are always implied. 
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criterion, determined by Monte Carlo simulation, and Nb is the number of back- 
ground pairs, determined from the data. This requires properly normalizing the 
number of signal pairs in the Monte Carlo to the number of signal pairs in the 
data. When the background becomes dominated by pairs with only one true D 
decay, we exclude from the background sample only those pairs in which both 
D candidates lie in a narrow range around the D mass. 

We iterate the procedure of finding the optimal selection criterion (always 
allowing variables to be reused in subsequent iterations) until the significance 
of the signal no longer increases. The selection criteria are optimized separately 
for each of five decay topologies of DD pairs: 2-2, 3-3, 2-3, 2-4 and 3-4, where 
each integer represents the number of charged particles in the decay. 5 We find 
that selection criteria are more often applied to one D candidate or the other, 
rather than to both, especially early in the optimization procedure. In several 
cases, a criterion will be applied to one of the D candidates, and a more stringent 
criterion involving the same discrimination variable will be applied to the other. 

After optimizing our selection criteria, we end up with a sample of 9254 events 
in the data with both D candidates in the mass range 1.7 to 2.0 GeV/c 2 and in 
the rapidity range —0.5 to 2.5. Only pairs in which the two D candidates have 
opposite charm quantum numbers are included in this sample. No significant 
signal for DD or D D pairs is observed. In Fig. 4, we plot the mass of the 
D candidate versus the mass of the D candidate, for all five types of pairs. 
Three types of candidate pairs are evident in this scatter plot. Combinatoric 
background pairs consisting of a fake D and a fake D candidate are spread over 
the entire plot. The density of these points decreases linearly with increasing 
candidate- D mass. Pairs containing one real and one fake D candidate appear 
as horizontal and vertical bands (called D and D ridge events, respectively). In 
the center of the plot, we see an enhancement due to the crossing of the two 
bands and due to real pairs of D mesons. 



5 Data Analysis 

In this section we describe the analysis procedures by which we determine the 
number of signal events in the full data sample shown in Fig. 4, as well as 
in each bin of the physics variables used to study the charm-pair production. 
Acceptance corrections include geometric acceptance, relative branching ratios, 
reconstruction efficiencies, and event selection efficiencies. 



5.1 Determination of Yields 

The experimental resolution for the D mass measurement in the E791 spectrom- 
eter depends on both the D decay mode and the xp of the D meson, and the 
mean reconstructed mass depends on the decay mode. Therefore, we fit to the 
normalized D mass defined as 

_ M - M D 

M n = , (6) 

cm 



5 We also searched for 4-4 pairs but the efficiency was too low to add much to the statistical 
significance of the sample. We did not use these pairs in the final analysis. 
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Fig. 4. Scatter plot of the D-candidate mass versus the D-candidate mass for 
the final unweighted charm-pair sample. 
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where M is the measured mass, Mp is the mean measured mass for the particular 
decay mode of the D candidate, and <tm is the experimental resolution for the 
particular decay mode and xf of the D candidate. 

In this analysis we use the maximum likelihood method which assumes we 
have N independent measurements of one or more quantities and that these 
quantities z are distributed according to some probability density function /(z; a.) 
where a is a set of unknown parameters to be determined. To determine the set 
of values a. that maximizes the joint probability for all events, we numerically 
solve the set of equations [40]: 

^|M = whcrc L (a) = n/C;a). 

The quantities that we measure for each event are the normalized mass of 
both the D and D candidate; i.e., z = (M,f ~" il7r , M,f +m27r ). The unknown 
parameters in the maximum likelihood fit are the number of signal events, Ns; 
combinatoric events, Nc] events with one real D and one combinatoric back- 
ground called D-ridge events, Nr>; events with one real D and one combinatoric 
background called Z?-ridge events, N-jj; the slope of the background K~m\'K dis- 
tribution, S D ; and the slope of the background if + TO 2 7r distribution, S D . That 
is, the unknown parameters are 

a = (N s ,N c ,N D ,N^,S D ,S B ). 

The terms K~m\Tt and K + m2^ refer to D or D decays into a kaon and mi 
pions. 

We construct our probability density function using the following two as- 
sumptions: (i) the normalized mass distribution for background K~m\ir and 
K + m 2 iT is linear in M% mi7T and M^ +m27! , and (ii) the normalized mass dis- 
tribution of real D's and real D's is Gaussian with mean of and sigma of 1. 
Under these assumptions, which are correct for our data, the probability den- 
sity functions — normalized to unity in the two-dimensional window defined by 
\Mn K mi7r )| < 6.5 and \Adi K m27T "> \ < 6.5 — for each class of events is 

Combinatoric background: P c = 1/169 + S D M*~ mi7T + S^il^f +m2Tr , 



D- Ridge background events: Pd = 
D- Ridge background events: Pjj = 
Signal events: Ps = 



(—k l NcLS D M K+m27T ) e~ {M * mi *)72 

[ 13V2^ + N- D 1V1 n ) e 

J_„-((M,f~ ml,r ) 2 + (M,f + m2 ' r ) 2 )/2 



The distribution for each set of events is NiPi. The overall probability density 
function is then 

. v = N C P C + N D P D + NjyPjy + N S P S 

S[Z,a) Nc + Nd+N-^ + Ns 
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In this analysis, we use the extended maximum likelihood method [41, pg. 
249] in which the number of DD candidates found, N D p, is considered to be 
one more measurement with a Gaussian probability distribution G(N D p, of 
mean fi = Nq + Nd + N-^ + N$ and a — ^/JL. Our likelihood function is then 

L = G(N d15 ,h) ]Jf(z t]a ). (7) 

i=i 

To maximize the likelihood, we use the function minimization and error anal- 
ysis FORTRAN package MINUIT [40]. Figure 5 shows the function 
N D j^f(z; a) that maximizes the likelihood function for the final sample of DD can- 
didates from Fig. 4 with \M n \ < 6.5, the mass range used for all fits in this 
analysis. The projections of the fit onto the D and D axes are compared to the 
data in Fig. 6. The projected background contains both ridge events (one real 
D and one combinatoric background) and events with two combinatoric back- 
ground candidates. Therefore, the background under the charm-pair signal in 
the projected distribution is a linear distribution plus a Gaussian distribution, 
shown as the dotted line in the figure. The net charm-pair signal is shown as the 
residual after background subtraction. 



5.2 Acceptance Corrections 

We determine the size of acceptance and smearing effects with a sample of ap- 
proximately 7000 Monte Carlo simulated pairs that pass the same selection cri- 
teria as the real data. The size of the resolution with which we measure each 
charm-pair physics variable is much smaller than the range over which we bin 
that variable. Therefore, we ignore smearing effects. 

We incorporate acceptance effects in the likelihood function for the fit by 
replacing the probability pi for event i by (pi) Wi where Wi is the weight for event 
i [42] . The weight is inversely proportional to the efficiency and is normalized 

N — — 

such that E i J^ D Wi — N D -^, where N D j=j is the number of DD candidates in the 
final sample. Corrections for relative branching fractions are also included in w% , 
as described below. By construction, the mean of the weights is equal to 1. The 
standard deviation of the weights is 1.3. The total number of DD events found 
in the unweighted fit is N s = 791 ± 44. For the weighted likelihood fit, we find 
N s = 910 ±45. 6 

The efficiency depends not only on the detector acceptance but also on the 
relative branching fractions for the detected decay modes. By correcting for 
branching fractions, the final efficiency-corrected distributions reflect the relative 
production rates of the four types of DD pairs (D°D°, D°D~, D + D°, and 
D + D~) rather than the relative detected rates. We use the values B(D + — > 

6 Since the sum of all weights is normalized to equal the number of DD candidates, the fact 
that the number of signal events is significantly larger for the weighted data sample indicates 
that, on average, the weights for signal events arc larger than for background events. Since 
we correct for relative efficiencies, not absolute efficiencies, the absolute number of weighted 
signal events has no significance. It is only of interest in interpreting the figures. 
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Fig. 5. The function that maximizes the joint probability of the unweighted 
charm-pair candidates from Fig. 4 with \M n \ < 6.5. The axes are the normalized 
candidate- D masses M n defined in Eq. 6 in the text. 
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Fig. 6. One-dimensional projections of the charm-pair normalized mass distribu- 
tions for D and D candidates, for unweighted events (top), and weighted events 
(bottom). The dashed curves are the fit projections. The dotted curves are the 
fit projections for the non-signal part of the fit. This includes background from 
one real D (D) and one fake D (D). The shaded histograms are the background- 
subtracted signals. 
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K-7T+7T+) = (9.1 ± 0.6)%, B(D° -» K-1T+1T-TT+) = (8.1 ± 0.5)%, and B(D° -» 
^-7r+) = (4.01±0.14)%[43]. 

A minimal independent set of properties that the acceptance could depend on 
is the decay mode of each of the D mesons (K~tt + , K~ir + ir + , or K ~7r + 7r~7r + ), 
the rapidity y, the transverse momentum p t , and the azimuthal angle <j> 01 eacn 
of the D mesons. In principle, we can use Monte Carlo simulated events to 
determine the acceptance for a particular candidate pair. The problem is the 
large number of Monte Carlo events that is needed to span such a large space (a 
48-dimensional space, six variables for each pair of decay modes used). However, 
the efficiency function can be factorized for each combination of the D decay 
modes, greatly increasing the statistical power of the Monte Carlo. 

Using the Monte Carlo simulated events, we find that the acceptance of the 
D is independent of the D, and vice versa, at the level of the statistical precision 
of the simulation. This is true in spite of the correlations in the selection criteria 
described in Section 4. For each one, however, the shape of the acceptance as 
a function of y depends on both the number of particles in the decay, no or 
n-jj, and, at high y, whether the candidate decay is a D or D. The shape of the 
acceptance as a function of pt depends only on the number of particles in the 
decay, or n-jy. It is also found that the acceptance does not depend on the 
azimuthal angle <j> of the D or D. Therefore, the acceptance function factorizes 
as follows: 

A(n D , rip, y D ,y^,p t ,D, P t ^, §d , <%) 

= b nD b n _c{y Dl y^)d nD {p^ D )d n _{p t - B ) 1 (8) 

where the subscripts, superscripts, and functional dependences of the terms b, 
c, and d are explicit, showing how the factorization is done. 

We next determine which of the variables that describe the candidate pair, 
and for which the acceptance is not uniform, are correlated in the originally 
generated Monte Carlo. Such correlations could affect the apparent acceptance 
from the Monte Carlo if we simply integrate over a variable that is correlated with 
the variable for which we are determining the acceptance. We find that the most 
significant correlations in the Monte Carlo generator are between the variables 
yu and yjy, where the acceptance is not uniform, and between the variables <j>D 
and <fi-Q, where the acceptance is uniform. Therefore, we cannot simply integrate 
over yo when determining the acceptance as a function of y-^. Instead, we use the 
Monte Carlo to determine the two-dimensional acceptance function c(yo, Ud) — 
which is independent of the Monte Carlo generated correlations — for each of the 
possible values of no and n-^ and use this function in Eq. 8. (This removes any 
dependence on the physics assumptions of the Monte Carlo from this equation.) 
Because of the uniform cj> acceptance for the observed events, the double variable 
technique is not needed for <j>. Finally, the weight Wi is calculated for each event 

i N — 

such that it is proportional to J.— and S- D , D Wi = N nT; , where N n jr is 

B{D)B(D)A uu uu 

the number of unweighted DD candidates in the final sample. Here B is B(D + — > 
K~tt + tt + ) for the charged D candidates and B is B(D° — ► K~ir + iT~ir + ) + 
B(D° — > K~ir + ) for the neutral D candidates except for D°D° events where 
B(D°)B(D°) = B(D° -» K-7T+) 2 +2B(D° -» K-ir+)B(D° -» K-ir+ir-ir+) 
due to the exclusion of 4-4 pairs from the final sample. 
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5.3 Checks &: Systematic Errors 

Sources of systematic errors in our measurements include effects associated with 
the fitting procedure used to obtain the yields, the finite statistics of the Monte 
Carlo data sample used to generate the acceptance contributions to the weights, 
and imperfections in our modeling of the apparatus in the Monte Carlo. 

For all the measured distributions, we compared the data to the two dimen- 
sional normalized mass distributions; in all cases, the fits qualitatively match 
the data. (For example, see Fig. 6.) 

We also checked the fitting procedure by comparing the yields with those 
given by a simple counting method. In this method the normalized mass scatter- 
plot was divided into regions corresponding to different combinations of signal, 
ridge, and combinatoric background events. The number of signal events was 
then found by subtracting the properly normalized number of events in the 
ridge and background regions from the central signal region. The results are in 
agreement, but the fitting technique gave smaller statistical errors, as expected. 

The effect of the finite statistics in the Monte Carlo was determined by 
repeating the fits for the yields in each kinematic bin while varying the weight 
of each event randomly according to a Gaussian whose width corresponded to 
the statistical error on the weight. The systematic errors on the yields generated 
by this process were about 20% of the statistical errors from the fit, which are 
negligible when added in quadrature. 

As demonstrated in Figs. 7 to 12, in most cases the weighted and unweighted 
distributions are very similar. Statistical errors associated with modeling the 
acceptance are most important for events with large weights, but the number 
of events with large weights is small. We checked the effect of large weights 
by generating distributions without the large weight events, with no significant 
change. The distributions were also examined with all Kititit candidates elimi- 
nated, the source of most of the events with a large acceptance correction. Again, 
the change did not significantly affect the distributions. 

Another potential source of systematic error is uncertainty in the modeling 
of the beam-induced inefficiency in the centers of the drift chambers. The inef- 
ficiencies increased as the run progressed, primarily affecting £>'s at large xf- 
Since most D's in the DD events are at low to modest xf, the drift chamber 
inefficiencies did not have a significant effect on the relative efficiencies used in 
this analysis. We checked this by comparing the experimental results that had 
been efficiency-corrected with Monte Carlo simulations corresponding to differ- 
ent parts of the run, and found no significant differences. 

In summary, systematic errors were found to be small relative to statistical 
errors. 

6 Results 

In this section, we present the background-subtracted, acceptance-corrected charm- 
pair distributions from the data and compare them to theoretical predictions. 
The Appendix contains an extensive discussion of the theoretical predictions for 
charm-pair distributions. As discussed in Sec. 4, all distributions — experimental 

http : //link. springer .de/link/ service/ journals/10105/ index.html 



EPJdiroct CI (2008) Springer- Verlag 



91 



5 
4 

a. 
x 

5 3 
Z 

5 2 

T— I 

1 





1 



"£10 



cio 



10 



0.2 


0.4 






• 




r •*• 








r ^ 


t 



2.5 5 7.5 10 
P? [GeV 2 ] 



1 r 

0.8 



N/dy 


0.6 






| 


0.4 








0.2 



5 r 

-e- 

Z 3-, 

T3 



90 180 270 360 
<|> [deg] 



Fig. 7. Single-charm distributions for xp, y, pi and <f>, obtained from summing 
the D and D distributions from our signal DD events; unweighted (o) and 
weighted (•) data are described in Sec. 5.2. 
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o Unweighted Data, • Weighted Data, Sircle — Charm Prediction 




Fig. 8. Charm-pair distributions for Axp = xf,d~x F -p and Sxp = xp^+x p -p. 

The uncorrected single-charm predictions ( ) for Axf and Sxf are defined 

in Eqs. 9 and 10. Unweighted (o) and weighted (•) data are described in Sec. 5.2. 



o Unweighted Data, • Weighted Data, Single — Charm Prediction 
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Fig. 9. Charm-pair distributions for Ay = Vd — Vb an( ^ = Vd + Vd- The 

uncorrelated single-charm predictions ( ) for Ay and Sy are defined in Eqs. 9 

and 10. Unweighted (o) and weighted (•) data are described in Sec. 5.2. 
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o Unweighted Data, • Weighted Data, Single-Charm Prediction 




Fig. 10. Charm-pair distributions for \Apf\ = \p^ D —p^—\ and Sp\ = p\ D +p^—. 

The uncorrelated single-charm predictions ( ) for Ap\ and Spl are defined 

in Eqs. 9 and 10. Unweighted (o) and weighted (•) data are described in Sec. 5.2. 
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Fig. 11. Charm-pair distributions for A<j> = (minimum of \4>d — <?%| and 360° — 
\<t>D — 4*1)1) an d £4> — {^d + 4>~d modulo 360°). The uncorrelated single-charm 

predictions ( ) for A(f> and E<f> are defined in Eqs. 9 and 10. Unweighted (o) 

and weighted (•) data are described in Sec. 5.2. 
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o Unweighted Data, • Weighted Data, Single-Charm Prediction 
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Fig. 12. Charm-pair distributions for p 2 D — = \pt,D + P t 7j| 2 an d M D p. The 
solid histogram shows the uncorrelated single-charm predictions. The dashed 
p 2 D — histogram assumes that xf,d and x p -^ are uncorrelated and that p 2 D and 

p 2 — are uncorrelated, but that 4>d and <^>yj are correlated as shown in Fig. 11. 
Unweighted (o) and weighted (•) data are described in Sec. 5.2. 



and theoretical — are obtained after excluding any events in which the center- 
of-mass rapidity of either D meson or either charm quark is less than —0.5 or 
greater than 2.5. 

For the experimental results, the acceptance-corrected distributions are ob- 
tained from maximizing the likelihood function with weighted events as discussed 
in Sec. 5.2. The uncorrected distributions are obtained from maximizing the un- 
weighted likelihood function; the total number of signal DD events found in the 
unweighted fit is N s = 791 ± 44. 

If the two charm mesons in each DD event are completely uncorrelated, then 
the charm-pair distributions contain no more information than the single-charm 
distributions. Before comparing the observed distributions to theoretical pre- 
dictions, we use two methods to determine whether there exist correlations in 
the data. In the first method, described in Section 6.1, we convolute acceptance- 
corrected single-charm distributions to predict what the charm-pair distributions 
would be if the D and D were uncorrelated. Comparing these single-charm pre- 
dictions, the measured charm-pair distributions provide one measure of the de- 
gree of correlation between the D and D. In the second method, described in Sec- 
tion 6.2, we look directly for correlations by examining several two-dimensional 
distributions. For example, by finding the number of signal DD events per yjj 
interval in several yjj intervals, we can determine whether the shape of the yp 
distribution depends on the value of yjy. In Section 6.3, we compare our exper- 
imental distributions to the theoretical predictions discussed in Sec. 2 and in 
the Appendix. In Section 6.4, we examine integrated production asymmetries 
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among the four types of DD pairs — D°D°, D°D~ , D + D° 1 and D + D~ — and 
compare our experimental results to the predictions from the Pythia / Jetset 
event generator. 

6.1 Single-Charm Predictions 

In Fig. 7 we show the measured single-charm distributions for xf, y, Vt and 
<j), as defined in Sec. 1. The single-charm distributions are obtained by fitting 
the two-dimensional normalized mass distributions for only those DD pairs in 
which the value of the single-charm variable for the candidate D (or D) is in 
the appropriate interval for each bin. In this way, the contribution to the single- 
charm signal from the D and D ridge events is excluded. The distributions 
shown in Fig. 7 correspond to single-charm mesons from DD pairs in which the 
center-of-mass rapidity of both charm mesons lies between —0.5 and +2.5. Each 
distribution shown in Fig. 7 is obtained by summing the D and D distributions. 
We have checked and found that the D and D distributions are the same within 
statistical errors. 7 

The vertical axis of each distribution gives the fraction of signal mesons per 
variable v interval, P(v) = j^^ff-, where the total number of signal D mesons 
No is simply twice the number of signal DD events. Only a very small fraction 
(0% - 3%) of the signal events lie outside any of the ranges used in Figs. 7-11. 

For each single-charm variable v = Xf, y, Pt, and we obtain two measured 
charm-pair distributions: the difference in v for the two -D's, Av = vd — vjj, and 
the sum of the v's for the two D's, Sv = vo + v-jy. ( Ac/) is defined to be the mini- 
mum of | <f>u — ?%| and 360° — \4>d — <?% | , and E(f> is defined to be 4>d + <?% modulo 
360°.) In Figs. 8-11, we compare these measured charm-pair distributions to the 
charm-pair distributions one would generate from the measured single-charm 
distributions assuming the D and D are completely uncorrelated, calculated as 
follows: 



where P(v) refers to the single-charm distributions shown in Fig. 7. 

This convolution cannot be done with previously reported inclusive single 
charm distributions [25] since the inclusive distributions contain events which 
are excluded from the charm-pair sample, for example, events in which the xf 
of the unobserved D is outside the acceptance or in which the unobserved charm 
particle is a charm baryon. 

If the D and D in the signal DD events are completely uncorrelated, then the 
measured charm-pair distributions for Av and Sv should agree with the single- 
charm predictions, because both the charm-pair and single-charm distributions 

7 This might not be the case if the experiment had greater statistics, or if the data sample 
extended to a higher region in xp . 





and 




S(Sv vjj v^)P(v D )P(v^)dv D dvp , 



(10) 
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are for D mesons with exactly the same restrictions on the rapidity of both D 
mesons in the event. With the exception of the A(f> distribution (Fig. 11), the 
measured distributions are quite similar to the uncorrelated single-charm pre- 
dictions, indicating both that the correlation between the D and D longitudinal 
momenta is small and that the correlation between the amplitudes of the D and 
D transverse momenta is small. The measured Axf and Ay distributions, how- 
ever, are somewhat more peaked near zero than the single-charm predictions, 
possibly indicating slight longitudinal correlations. 

Two other commonly used charm-pair variables are the square of the net 
transverse momentum of the charm pair, p 2 D — — \pt,D + P t r>\ 2 ' an d the in- 
variant mass of the charm pair, M D jy. The measured distributions and the un- 
correlated single-charm predictions for these two variables are shown in Fig. 12. 
Obtaining these single-charm predictions for p 2 t D — and M D p is slightly more 

involved than for the Av and Sv variables because p 2 D — and M D -^ are not 

linear functions of xf,d, x f~d> < hj^ Pt d> an< ^ p\ jy- Rather, in terms of these 
single-charm variables, 



pIdd = pId+pId + 2 \lp 2 t.DP 2 t ^ co <<pD - <%), and 



M dd = J 2M l + ZEdEjj - 2^p 2 D p 2 -cos( ( j )D - <%) - 



S X F ,D X F .D 



where the D meson energy E is y M|, + p\ + , and s is the square of the 
center-of-mass energy of the colliding hadrons. We obtain single-charm predic- 
tions by randomly generating 10 8 DD events in which all three variables (xf, 4> 1 
and p 2 ) for both D mesons from each DD event are selected independently and 
randomly from a probability density function that is flat within the bins shown 
in Fig. 7, and zero elsewhere. Each event is weighted by 

^P{xf,d)P{4>d)P{pI d )P{x f ^)P{<^)P{p^), 

where P(v) refers to the single-charm distributions shown in Fig. 7 and \J\ is the 
Jacobian determinant of the transformation from the complete and independent 
set of variables (xf,d, x fd> fin, <j>T>i Pt,D> an( ^P 2 j)) to the se t ( x f.d, x fd> 
^D' p2 t,DD> and M dd)- Specifically, \ J\ = 




l + t /¥co S (fe-^) 



^(t-^cos^-^)) 



The measured distribution for M fl ^ agrees quite well with the uncorrelated 
single-charm prediction. The measured distribution for p 2 — , however, is steeper 
than the uncorrelated single-charm prediction, indicating the presence of corre- 
lations between p ttD and p t jj- The dashed histogram in Fig. 12 demonstrates 



http : //link. springer . de/link/service/ journals/ 10105/ index.html 



EPJdiroct CI (2008) Springer- Verlag 



97 



that this lack of agreement is not due to the correlations between cj>D and (p-^ 
evident in the A<p distribution in Fig. 11. This latter prediction is obtained by 
assuming that x f .d and x F p are uncorrelated and that p 2 D and p 2 — are uncor- 
rected, but that (f>r> and are correlated as shown in Fig. 11. The correlations 
in p 2 — should reflect similar correlations in Ad>, p 2 n , and p 2 —. since p 2 — 

t.DD r ' y t,D' 1 t,D' t.DD 

is a function of these variables. The fact that the disagreement in Fig. 12 is 
not so readily explained is a sign that the correlations can be subtle, and that 
there are additional correlations among the variables. In the following section 
we investigate correlations between A<j>, Pt,D, and p f -^ in more detail. 

6.2 Two-Dimensional Distributions 

A direct method for investigating whether the variables vd and vjj are correlated 
is to determine the number of DD signal events per vd interval for a series of vjy 
intervals. Such two-dimensional distributions show whether the v D distribution 
depends on the value of v^, and vice-versa. Given the limited number of DD 
pairs, we use coarse binning to see statistically meaningful effects. In Figs. 13, 
14 and 15, we show the results for v = xf 1 y and pf, respectively. In part (a) of 
each figure, we show the number of acceptance-corrected DD signal events recon- 
structed in nine (vd, fjy) bins — three vd bins times three vjy bins. Note that 
the three bin sizes are not equal. From the information in this two-dimensional 
plot, several normalized one-dimensional plots are created, facilitating our abil- 
ity to detect differences in the shapes of the distributions. In particular, plot (b) 
in each figure shows the v D distribution for each vjy bin, g^j/^Vi, where N s is 
the number of events in the relevant bin, and iVj is the total number of events in 
the three vd bins. This normalization is chosen so that the integral over each vd 
distribution equals one. The symbols are defined in plot (a). Similarly, plot (c) 
in each figure shows the normalized vjj distribution for each vd bin. If the three 
sets of points in the figures (b) and (c) are statistically consistent, there are no 
significant correlations. Lastly, plots (d)-(f) simply rearrange the information 
shown in (b) and (c). Plot (d) shows the normalized vd and vjy distributions 
for the first vjj and vd bin, respectively; plot (e) shows results for the second 
bins; and plot (f) shows results for the third bins. In (d)-(f), agreement of the 
two sets of points implies that correlations in vjy are the same as correlations in 
vd ■ The two-dimensional plots show the actual number of acceptance-corrected 
DD signal events in each bin, whereas the one-dimensional plots, proportional 
to take into account the variation in bin size. 

Figure 13 indicates some correlation between Xf,d and x F p. In particular, 
the first-bin distributions are more peaked at low x p than the second- and third- 
bin distributions of figures (b) and (c). This result is consistent with Fig. 8, 
discussed above, which shows that the measured Axp distribution is somewhat 
steeper than the single-charm prediction. Because xp and y are highly correlated, 
Fig. 14 shows the same trends as Fig. 13. Figure 15 indicates that p 2 D and p 2 — 

are also slightly correlated — the second-bin p 2 D and p 2 — distributions are 
enhanced in the second bin; and the third-bin p\ D and p 2 — distributions are 
enhanced in the third bin of figures (b) and (c). This result should be compared 
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Fig. 13. (a) Number of weighted DD signal events dN s found in nine (xf,d, 
Xp-p) bins, (b) xf.d distribution (dN s /dxF,D)/Ni for each bin, where Ni 
is the total number of events in the three xf.d bins, (c) x F p distribution for 
each xf,d bin. (d)-(f) xf.d ( x fd) distribution for the first, second and third 
x f~d ( x f,d) bins, respectively. Open symbols show the xf,d distributions; closed 
symbols show the x p -^ distributions. The weighting procedure is described in 
Sec. 5.2. 
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Fig. 14. (a) Number of weighted DD signal events dN s found in nine (yo, 
y-p) bins, (b) yo distribution (d7V s /dyn)/A^ for each y-^ bin, where Ni is the 
total number of events in the three yr> bins, (c) y-^ distribution for each yu 
bin. (d)-(f) yo (j/jj) distribution for the first, second and third y-^ (yn) bins, 
respectively. Open symbols show the yo distributions; closed symbols show the 
y-p distributions. The weighting procedure is described in Sec. 5.2. 



with Fig. 10, which is consistent with no correlation. Correlations are also seen 
in Fig. 12 which shows that the measured D — distribution is somewhat steeper 

than the uncorrelated single-charm prediction. In all three figures (13-15), the 
shapes of the three vd distributions are remarkably similar to the shapes of the 
respective v-^ distributions as seen in figures (d)-(f). In Fig. 16 we investigate 
whether the separation in azimuthal angle between the D and D is correlated 
to the amplitude of the transverse momenta of the D and D. In particular, we 
determine the number of signal DD events per Atfi interval in Spj intervals and 
the number of signal DD events per A<fr interval in \A-p\ \ intervals. Although 
we find no significant correlation between A(f> and we find that A(f> and 

Spl are quite correlated. The A<j> distribution is more peaked at large Spl and 
the Spl distribution is flatter at large A(f>. A theoretical explanation for these 
correlations is discussed in the following section. 



6.3 Comparisons with Theory 

In this section, we compare all the acceptance-corrected distributions discussed 
in the previous two sections (Figs. 7-16) to three sets of theoretical predictions: 
the distributions of cc pairs from a next-to-leading-ordcr (NLO) perturbative 
QCD calculation by Mangano, Nason and Ridolfi (MNR)[3, 4]; the distributions 
of cc pairs from the Pythia/Jetset Monte Carlo event generator[5] which uses 
a parton-shower model to include higher-order perturbative effects [6]; and the 
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Fig. 15. (a) Number of weighted DD signal events dN s found in nine (p 2 iD , P 2 jj) 
bins, (b) p 2 D distribution (dN s /dp 2 D /Nj) for eachpj*— bin, where is the total 
number of events in the three D bins, (c) p^— distribution for each p\ D bin. 
(d)-(f) p\ D (p^-p) distribution for the first, second and third p 2 t — {p^ D ) bins, 
respectively. Open symbols show the p\ D distributions; closed symbols show the 
p 2 distributions. The weighting procedure is described in Sec. 5.2. 
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Fig. 16. (a) Number of weighted DD signal events diV s found in nine {A<j>, Spf) 
bins, (b) A(j> distribution (dN s /dA(f>)/Ni for each Ep\ bin, where Ni is the total 
number of events in the three A<f> bins, (c) Spf distribution for each A<f> bin. (d) 
Number of weighted DD signal events dA^ s found in nine (A<f>, \Ap1\) bins, (e) 
A(f> distribution for each \ ApJ\ bin. (f) \Ap\\ distribution for each A(f> bin. The 
weighting procedure is described in Sec. 5.2. 
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distributions of DD pairs from the Pythia (Version 5.7)/Jetset (Version 7.4) 
Monte Carlo event generator [5] which uses the Lund string model to transform 
cc pairs to DD pairs [7]. For all theoretical predictions, we use the default param- 
eters suggested by the respective authors, which are discussed in the Appendix. 
All distributions are obtained after excluding any candidates in which the center- 
of-mass rapidity of either the D or D (or, for the MNR calculation, the c or c) 
is less than —0.5 or greater than 2.5. 

6.3.1 Single-Charm Distributions Lack of agreement between an experi- 
mental charm-pair distribution and a theoretical prediction can arise if the the- 
ory does not model the correlations between the two charm particles correctly. 
However, it can also arise if the theory models the correlations correctly but 
does not correctly model the single-charm distributions. Hence, before compar- 
ing the experimental charm-pair distributions to theory, we first compare the 
acceptance-corrected single-charm distributions (xf, y, Vt an d 4>) to theory in 
Fig. 17. 

For the longitudinal momentum distributions — xf and y — the experi- 
mental results and theoretical predictions based on the default parameters do 
not agree. The experimental distributions arc most similar to the NLO and 
Pythia/ Jetset cc distributions, but are narrower than all three. The differ- 
ence between the Pythia/ Jetset cc and the Pythia/ Jetset DD longitudinal 
distributions shows the effect of the hadronization scheme that color-attaches 
one charm quark to the remnant beam and the other to the remnant target, 
broadening the longitudinal distributions. 

The experimental p% distribution agrees quite well with all three theoretical 
distributions. However, the Pythia/ Jetset cc distribution is somewhat flatter; 
the Pythia/ Jetset DD distribution is somewhat steeper. As expected, both 
the theoretical and experimental <j) distributions are consistent with being flat. 

6.3.2 Longitudinal Distributions for Pairs The experimental and theoret- 
ical Axf and Exf distributions (Fig. 18) and Ay and Ey distributions (Fig. 19) 
do not agree with theoretical predictions derived from default parameters. This 
may be a reflection of the disagreement between the measured single-charm lon- 
gitudinal distributions and theoretical models. As with the single-charm distri- 
butions, the experimental results are much closer to the two cc predictions than 
to the Pythia/ Jetset DD prediction, but narrower than all three predictions. 
The Pythia / Jetset hadronization scheme introduces a strong correlation be- 
tween the D and D which significantly broadens the Ay distribution. That is, 
hadronization tends to pull the DD apart, due to color string attachment to the 
incident hadronic remnants. As we show in Fig. 20, the Pythia/ Jetset DD 
Ay distribution is broader than the prediction we obtain by using the predicted 
single-charm distributions and assuming they are uncorrelated, as calculated us- 
ing Eqn. 9. In contrast, the experimental Ay distribution is slightly narrower 
than its uncorrelated single-charm prediction (Fig. 9). 

6.3.3 Transverse Distributions for Pairs In Figs. 21-23, we compare ex- 
perimental distributions to theoretical predictions for the following transverse 
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• Weighted Data, 



NLO cc, PYTH/JETcc, -- PYTH/JET DD 




2.5 5 7.5 10 

P? [GeV 2 ] 



90 180 270 360 
<t> [deg] 



Fig. 17. Single-charm distributions for xf, y, v\ an d <j>' weighted data (•) as de- 
scribed in Sec. 5.2; NLO QCD prediction ( ); Pythia/Jetset charm quark 

prediction ( ); and Pythia/ Jetset D meson prediction ( ). All distri- 
butions are obtained by summing the charm and anticharm distributions from 
charm-pair events. 
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Fig. 18. Charm-pair Axp and E x p distributions; weighted data (•) as described 
in Sec. 5.2; NLO QCD prediction ( ); Pythia/ Jetset charm quark predic- 
tion ( ); and Pythia/ Jetset D meson prediction ( ). 



• Weighted Data, NLO cc, - - - - PYTH/JET ad, - PYTH/JET DD 




Fig. 19. Charm-pair Ay and Ey distributions; weighted data (•) as described in 

Sec. 5.2; NLO QCD prediction ( ); Pythia/Jetset charm quark prediction 

( ); and Pythia/Jetset D meson prediction ( ). 
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Fig. 20. Pythia/Jetset DD prediction for Ay compared to the 
Pythia/Jetset single-charm prediction that is obtained by assuming 
that the D and D mesons are completely uncorrected (see Eq. 9). 



variables: |^lp 2 |, £p 2 , A<p, E<p, and p 2 D —. Any observed discrepancy between 
theory and data for the \Ap 2 \ , A<j>, and p 2 D — distributions is noteworthy because 

the single-charm pi and tfi experimental distributions agree well with theory. An 
observed discrepancy, therefore, must derive from the theory modeling the cor- 
relation between p t ,D and p t incorrectly. 

If pt.D arid p t -p were completely uncorrelated, then the single-charm predic- 
tions (Figs. 10-12) would provide good estimates for these three distributions. 
At the opposite extreme, if Pt,D and p t were completely anticorrelated — as 
in the leading-order perturbative QCD prediction — then the Ap\ distribution 
would be a delta function at An? = GeV 2 ; the p 2 „— distribution a delta func- 

1 t,DD 

tion at p 2 D — — GeV 2 ; and the A<f> distribution a delta function at A<p ~ 180°. 
Both the experimental distributions and the three sets of theoretical predictions 
lie between these extremes. None of the three experimental distributions, how- 
ever, is as steep as any of the theoretical predictions. The next-to-leading-order 
predictions are the steepest — that is, the next-to-leading-ordcr calculation pre- 
dicts the most correlation between p tjC and p t c as a model for p t ,D and p t 
Thus hadronization and higher-order perturbative effects smear out the cc cor- 
relations. The Pythia/Jetset hadronization scheme broadens the A<p distri- 
bution, bringing it closer to the experimental results. The same hadronization 

scheme also narrows the p 2 — and Api distributions, which makes them disagree 
^t.DD ^ l ' to 

even more with the experimental results. One mechanism which would broaden 
the p 2 D — and Ap\ distributions as well as the A(f> distribution (bringing all 

into better agreement with the experimental results) is to increase the intrinsic 
transverse momentum of the colliding partons in the beam and target hadrons 
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• Weighted Data, NLO cc, PYTH/JET cc, -- PYTH/JET DD 




2 4 6 2.5 5 7.5 10 



|Ap 2 | [GeV 2 ] Ep 2 [GeV 2 ] 

Fig. 21. Charm-pair |Z\pj | and Spf distributions; weighted data (•) as described 
in Sec. 5.2; NLO QCD prediction ( ); Pythia/ Jetset charm quark predic- 
tion ( ); and Pythia/ Jetset D meson prediction ( ). 

(see the Appendix, Sees. A. 5 and A. 6). An improved theoretical understanding 
may involve adding terms higher than NLO to calculations, although other au- 
thors find good agreement by choosing appropriate values for nonperturbative 
parameters [44]. 

6.3.4 Charm-Pair Invariant Mass In Fig. 23, we compare the experimental 
charm-pair invariant mass distribution to the Pythia/ Jetset prediction. The 
experimental distribution is steeper than the theoretical predictions. This 
is similar to the experimental single-charm xp (or y) distributions, which are 
also steeper than the theoretical predictions. (See Fig. 17.) 

In addition, the correlations introduced by the Pythia/ Jetset hadroniza- 
tion scheme broaden the invariant mass distribution. 

6.3.5 Two-Dimensional Distributions In Figs. 24-27, we examine the 
same two-dimensional experimental distributions discussed in Section 6.2. We 
now compare these experimental results to the three sets of theoretical predic- 
tions. In each figure, the top row shows the NLO perturbative QCD cc predic- 
tion; the middle row, the Pythia/ Jetset cc prediction; and the bottom row, 
the Pythia/ Jetset DD prediction. The experimental data points and errors 
are repeated in each row. 

The longitudinal distributions, xp and y, are shown in Figs. 24 and 25. 
The three theoretical predictions are quite different. The NLO cc predictions 
show no significant correlation between xf,d and x F -p (or between yjy and j/jj) 
and the xp,D and x p -^ distributions are quite similar. The Pythia/ Jetset 
cc predictions show a slight correlation and the xf.d and x p -^ distributions 
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• Weighted Data, 
30 r 



NLO cc 




45 90 135 180 90 180 270 360 
A<|> [deg] £<j) [deg] 



Fig. 22. Charm-pair A<j> and £<j> distributions; weighted data (•) as described in 

Sec. 5.2; NLO QCD prediction ( ); Pythia/Jetset charm quark prediction 

( ); and Pythia/Jetset D meson prediction ( ). 



• Weighted Data, NLO cc, PYTH/JETcc, -- PYTH/JET DD 



1.25 




2.5 5 7.5 10 4 6 



Pt 2 DD [GeV 2 ] M d5 [GeV] 

Fig. 23. Charm-pair D — and M D -^ distributions: weighted data (•) as de- 
scribed in Sec. 5.2; NLO QCD prediction ( ); Pythia/Jetset charm quark 

prediction ( ); and Pythia/ Jetset D meson prediction ( ). 
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are somewhat different. Due to the Pythia/ Jetset hadronization scheme, the 
Pythia/Jetset DD prediction shows the strongest correlation between Xf,d 
and Xp-p. Interestingly, in the Pythia/ Jetset DD prediction, Xf,d and x p -^ 
are anticorrelatcd; in the Pythia/Jetset cc prediction they are positively cor- 
related. The correlation patterns in the experimental results, although inconsis- 
tent with any of the theoretical predictions, are closest to the Pythia/Jetset 
cc predictions. 

In Fig. 26, we investigate the correlations between p\ D and P^-p- The three 
theoretical predictions and the experimental results all show similar trends. Al- 
though all the distributions are broader than the leading-order perturbative 
QCD prediction — a delta function at p\ D =p l — — they all shows signs of an 

enhancement in the p\d = p 2 -^ bins. The Pythia/Jetset cc distributions and 

the Pythia/Jetset DD distributions are very similar and resemble the exper- 
imental results more so than the NLO cc distributions. All of the theoretical 
third-bin distributions are significantly flatter than the experimental third-bin 
distributions. In contrast to the longitudinal distributions, all the pf D distribu- 
tions are very similar to the respective p"^— distributions. 

In Fig. 27, we investigate correlations between A<j> and Spt- For the A<j> de- 
pendence, a leading-order perturbative QCD calculation predicts a delta function 
at Acj) = 180°. We expect perturbative predictions to be more accurate as the 
energy scale Q of the partonic hard scattering increases: 



That is, we expect the A<j> distribution to be more peaked at 180° for DD events 
with larger Sp\. This behavior is clearly evident in our experimental distribu- 
tions as well as in all three theoretical predictions. The theoretical A<f> distribu- 
tions, however, for all three ZJpf bins, are significantly steeper than the respective 
experimental distributions. The NLO cc A<j> distributions are the steepest. The 
experimental and theoretical Spl distributions are in fairly good agreement, 
with the Spl distribution broadening as Acj) increases. No significant correlation 
between A<f> and \Ap%\ or between A<p and Ay is seen in the data or theory. 

6.4 Dependence of Yields and Longitudinal Correlations on Type 
of DD Pair 

6.4.1 Relative Yields In Table 3, we compare the experimental yields for 
each type of DD pair to the predictions from the Pythia / Jetset event gener- 
ator and to a naive spin-counting model. The experimental results are obtained 
by maximizing the weighted likelihood function where the weights account for 
both acceptance effects and the relative branching fractions of the reconstructed 
decay modes (Sec. 5.2). Again, for both data and Pythia/Jetset predictions, 
the results are for pairs in which the rapidities of both the D and D lie between 
—0.5 and 2.5. The Pythia/Jetset and naive spin-counting models both as- 
sume that vector D* production is three times more likely than pseudoscalar D 




(11) 
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Fig. 24. (a) The NLO QCD prediction for the number of cc events in nine (xf,d, 
x F jj) bins, normalized such that the number of NLO events equals the number 

of weighted DD signal events, (b) Experimental xf,d distribution for each x F p 
bin compared to the NLO QCD predictions. Each Xf,d distribution is normal- 
ized such that the integral over xf,d equals one. (c) Same as (b) for the x F -p 
distributions, (d)-(f) Same as (a)-(c) for the Pythia/ Jetset cc prediction, (g)- 
(i) Same as (a)-(c) for the Pythia/ Jetset DD prediction. Symbols represent 
weighted data; histograms represent theoretical predictions. A and — — corre- 
spond to the low bin; • and to the middle bin; □ and to the high 

bin. 
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Fig. 25. (a) The NLO QCD prediction for the number of cc events in nine (yD, 
y-jj) bins, normalized such that the number of NLO events equals the number 
of weighted DD signal events, (b) Experimental yo distribution for each y-^ bin 
compared to the NLO QCD predictions. Each yp distribution is normalized such 
that the integral over yo equals one. (c) Same as (b) for the y-^ distributions. 
(d)-(f) Same as (a)-(c) for the Pythia/Jetset cc prediction, (g)-(i) Same as (a)- 
(c) for the Pythia/Jetset DD prediction. Symbols represent weighted data; 
histograms represent theoretical predictions. A and — — correspond to the low 
bin; • and to the middle bin; □ and to the high bin. 
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Fig. 26. (a) The NLO QCD prediction for the number of cc events in nine (j$ D , 
p 2 t — ) bins, normalized such that the number of NLO events equals the number 

of weighted DD signal events, (b) Experimental pf D distribution for each p 2 — 
bin compared to the NLO QCD predictions. Each p\ D distribution is normalized 
such that the integral over p\ D equals one. (c) Same as (b) for the p 2 — distribu- 
tions. (d)-(f) Same as (a)-(c) for the Pythia/ Jetset cc prediction, (g)-(i) Same 
as (a)-(c) for the Pythia/Jetset DD prediction. Symbols represent weighted 
data; histograms represent theoretical predictions. A and - - correspond to 
the low bin; • and to the middle bin; □ and to the high bin. 
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Fig. 27. (a) The NLO QCD prediction for the number of cc events in nine (A<fi, 
Spf) bins, normalized such that the number of NLO events equals the number of 
weighted DD signal events, (b) Experimental A<j> distribution for each Spl bin 
compared to the NLO QCD predictions. Each A<p distribution is normalized such 
that the integral over A<j> equals one. (c) Same as (b) for the Epf distributions. 
(d)-(f) Same as (a)-(c) for the Pythia/ Jetset cc prediction, (g)-(i) Same as (a)- 
(c) for the Pythia/ Jetset DD prediction. Symbols represent weighted data; 
histograms represent theoretical predictions. A and — — correspond to the low 
bin; • and to the middle bin; □ and to the high bin. 
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production due to the number of spin states and that contributions from higher 
spin states are negligible. They also use the known D r± branching fractions, 
B(D*+ -> D°tt+) = 68.3% and B(D*+ -» D+X) = 31.7%, to determine D pro- 
duction. The differences between the Pythia/Jetset and naive spin-counting 
model come from the Pythia/Jetset hadronization scheme — in particular, 
the rate of coalescence. As discussed in Sec. 2.2, a charm quark tied to a valence 
quark by a low-mass string can coalesce with that valence quark into a meson. 
This will tend to increase the rate of D~(cd) and D°(cu) production in the for- 
ward region for the E791 n~(ud) beam. Since D*~ decays to D°, production 
of D° will also be enhanced. This effect is seen in Tabic 3 where the number 
of pairs that contain a D + is reduced while the number of pairs that contain a 
D~,D°, or D° is increased in the Pythia/Jetset model relative to the naive 
model. Both models agree with data as far as the relative ordering but predict 
too many D°D° pairs and too few D + D~ pairs. 



Table 3. Normalized acceptance-corrected experimental yields for the four types 
of DD pairs compared to predictions from a naive spin-counting model and the 
Pythia/Jetset event generator. Experimental and Pythia/Jetset yields are 
for —0.5 < y D j) < 2.5. Statistical and systematic errors are given for the data. 





Data 


Spin-counting model 


Pythia/Jetset 




0.50 ±0.04 ±0.01 


0.572 


0.604 ±0.002 


D°D~ 


0.20 ±0.02 ±0.01 


0.184 


0.208 ±0.002 


D+D° 


0.18 ±0.02 ±0.01 


0.184 


0.138 ±0.002 


D+D- 


0.12 ±0.02 ±0.01 


0.060 


0.051 ±0.001 



6.4.2 Correlations Between the D and D Longitudinal Momenta As 

shown in Fig. 28, in the Pythia/Jetset hadronization scheme, the correlation 
between yp and yjj is quite different for each of the four types of DD pairs. 
In Fig. 29, we investigate whether this is also true for data. Given the limited 
size of our data sample, we can only search for gross asymmetries in the (j/d, 
y-jj) distributions. We obtain the four plots in Fig. 29 by bisecting the two- 
dimensional (yo, y-jj) distributions along the following four lines (v = a): Ay = 0, 
Sy = 1.2, yn = 0.6, and y-^ — 0.6. These four lines are indicated by dashed 
lines in Fig. 28. 

To search for possible differences in asymmetries, we determine whether the 
fraction of signal events on one side of a given line depends on the type of DD 
pair. Specifically, for both theory and data, we show in Fig. 29 

_ N i {v>a)-N i {v<a) 
AvW ~ Ni(v > a) + Ni(v < a) [ ' 

where i = {D°D ,D a D-,D + D°,D + D-) and JV; is the number of signal DD 
events of type i. The Pythia/Jetset A Sy distribution is fairly flat, indicat- 
ing no significant differences among the four DD types for the Sy distribu- 
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Fig. 28. Pythia/Jetset prediction for the (y D , y-p) distribution for each of the 
four types of DD pairs. The dashed lines help define the asymmetry functions 
A v {i) (Eq. 12) shown in Fig. 29. 
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Fig. 29. The asymmetry function A v (i) = ^^>°)^^<°j , where i =(D°D°, 
D°D~, D + D°, D + D~) and Ni is the number of signal DD events of type i, 
for (v, a) = (Ay, 0), (Ey, 1.2), (j/d, 0.6), and (y-p, 0.6). Both the weighted 
data (open symbols) and the Pythia/Jetset (closed symbols) distributions 
correspond to events in which —0.5 < yn, < 2.5. 
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tion. The Pythia/ Jetset A^ y , A Vdi and A y _ distributions, however, indi- 
cate significant differences, all of which are easily interpreted in terms of the 
Pythia/ Jetset coalescence mechanism discussed in Sec. 2. Unfortunately, our 
experimental errors are of the same order as the degree of differences in the 
Pythia/ Jetset predictions. The experimental A VD distribution, for example, 
is consistent with the Pythia/ Jetset prediction, but it is also consistent with 
being flat. Similarly, the experimental A^ y distribution is fairly consistent with 
the flat Pythia/ Jetset prediction, but it also shows some indication of a dif- 
ference between D°D° and D + D~ . The most significant difference between the 
experimental results and the Pythia/ Jetset predictions occurs for the DD 
types D°D° and D + D~ in the Ay_ distribution. Both theory and data indi- 
cate a difference between D°D° and D + D~: however, we find experimentally 
that Ay_(D°D°) > A y _(D + D~), whereas the Pythia/ Jetset model finds 
A yj _(D+D-) > Ay_(D°W). 

7 Conclusions 

We fully reconstructed 791 ± 44 true DD pairs after all background subtrac- 
tions. This is the largest such sample of charm pairs used in an analysis of the 
hadroproduction of cc to date. The full reconstruction of the final states of both 
D mesons offers several advantages over some of the previous studies that have 
used partially reconstructed D candidates. We are able to correct the data for 
both detector inefficiencies and for the branching fractions of the observed decays 
so that the acceptance-corrected distributions represent the produced mixture 
of D mesons, rather than the detected mixture. Because the final states are 
fully reconstructed, we are able to calculate both the magnitude and direction 
of the D momenta. Therefore, we are able to thoroughly investigate the degree 
of correlation between both the transverse and longitudinal components of the 
momenta, with respect to the beam direction, of the D and D. 

We have compared all the measured acceptance-corrected distributions to 
predictions of the fully differential next-to-leading-order calculation for cc pro- 
duction by Mangano, Nason and Ridolfi [3, 4], as well as to predictions from 
the Pythia / Jetset Monte Carlo event generator [5] for cc [6] and DD produc- 
tion [7]. 

7.1 Transverse Correlations 

Our measurements indicate that the transverse momenta of the D and D in 
charm-pair events are correlated in several ways. (See Sees. 6.2 and 6.3.5.) The 
square of the amplitudes of the D and D transverse momenta are slightly cor- 
related (Fig. 15). The directions of the D and D in the plane transverse to the 
beam axis are significantly correlated (Fig. 11). The separation in azimuthal 
angle, A(f>, is significantly correlated with the sum of the squares of the D and 
D transverse momenta, Spl (Fig. 16). These features have been observed by 
several other experiments [8, 9, 10, 11, 12, 13, 14, 15, 16]. Using the default 
parameters, the three models yield the same trend in correlations as we find in 
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data. The models also predict that the relative transverse angles of the D and 
D are more correlated than we find in data (Fig. 22). These results provide an 
opportunity to tune the default parameters, or add higher order terms in the 
models, to obtain better agreement with data [44]. 

7.2 Longitudinal Correlations 

Our measurements indicate that the longitudinal momenta of the D and D from 
charm-pair events are slightly correlated. The measured Axf and Ay distri- 
butions (Figs. 8-9) are somewhat narrower than what one would predict from 
the observed single-charm predictions assuming no correlations. The xf,d (ud) 
distribution depends on the value of x F jj (yjj), and vice- versa (Figs. 13-14). 

The single-charm x p and y distributions from the three theoretical models do 
not agree with each other or with the measured distributions (Fig. 17). The three 
models predict different correlations between the charm and anticharm longitu- 
dinal momenta (Figs. 24-25) — the next-to-leading-order calculation predicts 
no significant correlation; the Pythia/ Jetset cc prediction indicates a slight 
positive correlation; and the Pythia/ Jetset DD prediction indicates a strong 
negative correlation. The DD data agree best with the Pythia / Jetset cc pre- 
diction. The disagreement between the models and the data might be corrected 
by adjusting the non-perturbative parameters in the models, or by adding higher 
order terms. 

7.3 Dependence of Yields and Longitudinal Correlations on Type 
of DD Pair 

The relative yields of the four types of charm pairs, D°D°, D°D~ , D + D°, and 
D + D~, as calculated in the Pythia/ Jetset event generator, agree with data 
as far as their ordering but disagree with regard to number of pairs produced, 
predicting too many D°D° pairs and too few D + D~ pairs. (See Table 3.) We 
studied the degree to which longitudinal correlations depend on the type of 
DD pair in data and in the Pythia/ Jetset event generator. Although we 
see differences between data and the event generator (Fig. 29), the statistical 
uncertainties on the measured correlations are too large to make any conclusive 
statements. 

7.4 Summary and Discussion 

The charm-pair distributions presented in this paper provide an opportunity 
to extend our understanding of charm production beyond what was previously 
possible with single-charm and lower-statistics or partially reconstructed charm- 
pair distributions. The measured distributions and observed correlations can be 
compared to the predictions of models, testing assumptions in the models and 
providing discrimination among different values for the free parameters in the 
models. Some comparisons have been made in the paper and in the Appendix. 

Before comparing the measurements to predictions, we looked for correlations 
directly in the charm-pair data by comparing the observed DD pair distributions 

http : //link. springer . de/link/service/ journals/ 10105/ index.html 



EPJdiroct CI (2008) Springer- Verlag 



118 



with the convolution of the measured single-charm distributions assuming no 

correlations. We find that the charm-pair distributions are quite similar to the 

convoluted single-charm distributions, indicating little correlation between the 

two charm mesons in an event, with the exception of the distributions for AS 

and p 2 —. The AS distribution shows clear evidence of correlations, and the 
t.DD r ' 

p D — distribution is steeper than the uncorrelated single-charm prediction. In 
addition, the data are consistent with possible small correlations in the Axp 
and Ay pair distributions, which are somewhat more peaked near zero than the 
single-charm convolutions. 

In the comparisons of the measured and predicted charm-pair distributions, 
we observe less correlation between transverse momenta and different correla- 
tions between longitudinal momenta than theoretical models predict, for the 
default values of parameters in the models. Work by other authors suggests a 
different set of parameters might provide better agreement [44] . Both the single- 
charm and charm-pair distributions agree best with the predictions for charm 
quark (rather than D meson) production, possibly caused by an accidental can- 
cellation of color-dragging and fragmentation effects. Also, the A6 distribution 
is more similar to the prediction of the NLO theory at higher Up 2 . 

In the Appendix, we investigate the sensitivity of single-charm and charm- 
pair distributions to various theoretical assumptions. We conclude that the pre- 
dictions depend not only on unknown parameters such as the mass of the charm 
quark and the intrinsic transverse momentum of the partons that collide to form 
the cc pair, but also on the values of the renormalization and factorization scales. 
The measurements reported here, and the charm-pair measurements from pho- 
toproduction experiments, should allow the free parameters in the theoretical 
models to be further constrained. 
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A Theoretical Predictions for Charm-Pair Distributions 

In Sec. 2, we introduced the theoretical framework used to describe the hadropro- 
duction of DD pairs. In particular, we discussed: 

• the leading-order perturbative QCD description of the hadroproduction of 
cc pairs; 
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• higher-order perturbative corrections to the leading-order calculation; 

• the addition of intrinsic transverse momentum to the hard-scattering par- 
tons that collide to form the cc pair; and, lastly, 

• the hadronization of cc pairs to observable DD pairs. 

Using this framework, we investigate how sensitive single-charm and charm-pair 
distributions are to various theoretical assumptions. All predictions discussed in 
this Appendix are for a 500 GcV/c tt~ beam incident on a nuclear target — the 
same beam-target as the data from experiment E791 reported in this paper. 

The Pythia/ Jetset event generator depends on many parameters. Un- 
less otherwise mentioned, we use the default settings for all parameters. The 
next-to-lcading order perturbative QCD calculation depends on the following 
six parameters: 

• the mass of the charm quark, m c , 

• the beam and target parton distribution functions, f n and f N , respectively, 

• Aqcd, the free parameter that must be determined experimentally, which 
roughly defines the mass scale below which quarks and gluons do not be- 
have as independent, free partons — that is, below which perturbative 
QCD calculations are no longer valid, and 

• the renormalization and factorization scales, fin and Hf- 

The pairs of pion and nucleon parton distribution functions considered in this 
section, obtained from the CERN computer library package PDFLIB [45], are 
listed in Table 4. Parton distribution functions depend on the fraction x of the 
hadron momentum carried by the hard-scattering parton, on both the factor- 
ization and renormalization scales, and on Aqcd- In the parton distribution 
functions accessible from PDFLIB, the renormalization scale is defined to be the 
same as the factorization scale. 

For each parton distribution function listed in Table 4, we specify the square 
of the minimum factorization scale allowed, Hq\ whether the evolution equations 
were calculated to leading-order (LO) or to next-to-leading order (NLO) and 
the value of Aq CD used in the fit. Querying PDFLIB for the value of a parton 
distribution function at a scale below /jq gives undefined results. The default set 
of parton distribution functions for the Pythia/ Jetset event generator is set 
(5) in Table 4; the default suggested by the authors of HVQMNR is set (1). 

When possible, we choose pion and nucleon distribution functions that are fit 
assuming similar values for Aqcd- For all predictions shown below, the Aqcd 
used in the next-to-leading order calculation of the partonic cross section is de- 
fined to be the same as the Aqcd used to extract the nucleon parton distribution 
function f N . 

The degree to which the charm-pair distributions are sensitive to variations 
in and hf gives an indication of how important higher-order corrections are; 
that is, an indication of how much (or little) we can trust the calculation. 
In general, one tries to minimize higher-order contributions by choosing /ir and 
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Table 4. The pairs of pion and nucleon parton distribution functions considered 
in this Appendix, obtained from the CERN FORTRAN package PDFLIB. The 
functions have been extracted from fits to data assuming a fixed value of Aqc d ■ 
The functions are undefined below the minimum scale /i 



Set 




Name 


Mo [GeV 2 ] 


[MeV] 


Order 


Ref. 


(1) 


r 


SMRS-P2 


5 


190 


NLO 


[46] 




f N 


HMRS-B (4.90) 


5 


190 


NLO 


[47] 


(2) 




GRV-P 


0.3 


200 


NLO 


[48] 




jN 


GRV 


0.3 


200 


NLO 


[49] 


(3) 


r 


SMRS-P2 


5 


190 


NLO 


[46] 




f N 


HMRS-B (8.90) 


5 


100 


NLO 


[50] 


(4) 


r 


ABFKW-P3 


2 


281 


NLO 


[51] 




f N 


HMRS-B (8.90) 


5 


300 


NLO 


[50] 


(5) 9 


r 


OW-P1 


4 


200 


LO 


[52] 




f N 


CTEQ 2L 


4 


190 


LO 


[53] 



[if to be of the same order as the energy scale Q of the hard-scattering process. 
However, this scale cannot be defined unambiguously. One reasonable choice is 

q - yU+^fc d3) 

The default setting for the Pythia/ Jetset event generator is fiR = \lf = Q, 
leading to factorization scales as low as the mass of the charm quark, m c , which 
by default is set to 1.35 GeV. The parton distribution functions used by the 
Pythia/ Jetset event generator, however, are only defined for scales above 2 
GeV. This problem is handled by setting the parton distribution function to 
f(x, no) f° r au factorization scales less than ^o- 

The suggested default for the HVQMNR program is hr = Q an d Hf = 2Q. 
Given their suggested default for the mass of the charm quark of m c = 1.5 GeV, 
this choice ensures that the factorization scale will never go below the minimum 
allowed scale, = V5 GeV. 

In Figs. 30-36 we show single-charm and charm-pair distributions for a wide 
range of theoretical assumptions. When obtaining these theoretical predictions, 
we only allow charm-pair events in which both charm rapidities are greater than 
—0.5 and less than 2.5. For the HVQMNR generator, which does not hadronize 
the cc pair to charmed mesons, the cut is on the charm quark rapidities. For the 
Pythia/ Jetset generator, the cut is on the D meson rapidities. In Table 5, 
we show which generator (HVQMNR or Pythia / Jetset) and what theoretical 
assumptions are used in each figure. 
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The same set of single-charm and charm-pair distributions are shown in each 
figure. Each charm particle in a charm-pair event can be described using three 
variables. A common choice of independent variables for single-charm analyses 
is xf, v\i an d <fi- 

We ignore the latter variable because all theoretical predictions give a flat 
4> distribution. Although xp and y are very correlated, we show predictions for 
both distributions. For each single-charm variable v, we obtain predictions for 
two charm-pair distributions: Av = v c — Vc and Ev — v c + v-. (Atfi is defined to 
be the minimum of \<p c — <f>c\ and 360° — \<p c — <fe|.) As with the single-charm 4> 
variable, we ignore the charm-pair £<f> variable because all theoretical predictions 
give a flat £</> distribution. We do not, however, ignore the A(f> distribution which 
is very sensitive to theoretical assumptions. Two other commonly used charm- 
pair distributions that we examine are the square of the transverse momentum 
of the charm-pair, p 2 c - — \pt,c + Pt,c\ 2 , an d the invariant mass of the charm-pair, 
M- c . 

The vertical axis of each distribution shown in Figs. 30-36 gives the fraction 
of single-charm (charm-pair) events per variable v interval, jj^-, where N is 
the total number of single-charm (charm-pair) events generated. The number of 
single-charm events generated is, of course, just twice the number of charm-pair 
events generated. 

A.l Sensitivity to Higher-Order Perturbative Corrections 

In Fig. 30, we compare the complementary methods used by the HVQMNR pro- 
gram and the Pythia/Jetset event generator to include higher-order pertur- 
bative corrections to the leading-order partonic cross section. As discussed in the 
previous section, the Pythia/ Jetset event generator, beginning with leading- 
order matrix elements, uses parton showers to include higher-order perturbative 
effects, whereas the HVQMNR program calculates the next-to-leading order cc 
cross section. To more directly compare these two approaches, we change three of 
the default Pythia/Jetset settings — m c , f* and f N — to match the default 
HVQMNR settings. (See Table 5.) We obtain Pythia/Jetset cc distributions 
assuming no intrinsic transverse momentum for the interacting partons, as well 
as assuming o kt — 0.44 GeV which is the Pythia/Jetset default. As argued by 
T. Sjostrand, the intrinsic transverse momentum may, at least in part, be seen 
as a replacement for gluon emission that is truncated in the parton shower ap- 
proach due to the introduction of an energy scale below which the parton shower 
evolution is stoppcd[54]. In Fig. 30, we also show the HVQMNR leading-order 
distributions to emphasize which distributions are, and which are not, sensitive 
to higher-order corrections. 

Figure 30 shows that higher-order perturbative corrections do not signifi- 
cantly affect the shapes of most of the single-charm and charm-pair distributions. 
That is, the HVQMNR leading-order and next-to-leading predictions for all dis- 
tributions are very similar — except for the \ Ap\ |, A<j), and p 2 c - distributions. 
In the leading-order calculation, these latter distributions are delta functions — 
at (GeV/c) 2 , 180°, and (GeV/c) 2 , respectively — because the leading-order 
charm and anticharm quark are back-to-back in the plane transverse to the beam 
axis. 
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Table 5. The settings used by the HVQMNR and Pythia/Jetset generators 
to obtain the single-charm and charm-pair distributions shown in Figs. 30-36. 
The set of pion and nucleon parton distribution functions (PDF), labeled (1) 
through (5), arc defined in Table 4. A "Y" indicates that parton showers (PS) 
are included in the Pythia/Jetset event generator; an "N" indicates that they 
are not included. The last column describes the histogram style corresponding 
to the settings in that row, in the figure listed in the second to last column. 



Generator 




fur 




Hf/Q 


m c 
[GeVJ 


[GeVJ 


Fig./PS 




MNR 


NLO a 


(1) 


1.0 


2.0 


1.5 





30 


solid 


MNR 


LO 6 


(1) 


1.0 


2.0 


1.5 







dashed 


P/J 


cc c 


(1) 


1.0 


1.0 


1.5 





Y 


dotted 


P/J 


cc 


(1) 


1.0 


1.0 


1.5 


0.44 


Y 


solid 


MNR 


NLO 


(1) 


1.0 


2.0 


1.5 





31 


solid 


MNR 


NLO 


(1) 


1.0 


2.0 


1.2 







dashed 


MNR 


NLO 


(1) 


1.0 


2.0 


1.8 







dotted 


MNR 


NLO 


(1) 


1.0 


2.0 


1.5 





32 


solid 


MNR 


NLO 


(2) 


1.0 


2.0 


1.5 







dashed 


MNR 


NLO 


(3) 


1.0 


2.0 


1.5 







dotted 


MNR 


NLO 


(4) 


1.0 


2.0 


1.5 







solid 


MNR 


NLO 


(2) 


0.5 


0.5 


1.5 





34 


solid 


MNR 


NLO 


(2) 


1.0 


1.0 


1.5 







dashed 


MNR 


NLO 


(2) 


1.5 


1.5 


1.5 







dotted 


P/J 


DD d 


(5) 


1.0 


1.0 


1.35 


0.44 


35/ Y 


solid 


P/J 


cc 


(5) 


1.0 


1.0 


1.35 





Y 


dashed 


P/J 


cc 


(5) 


1.0 


1.0 


1.35 


0.44 


N 


dotted 


P/J 


DD 


(5) 


1.0 


1.0 


1.35 





N 


solid 


P/J 


DD 


(5) 


1.0 


1.0 


1.35 


0.44 


36/ Y 


solid 


P/J 


DD 


(5) 


1.0 


1.0 


1.35 


0.7 


Y 


dashed 


P/J 


DD 


(5) 


1.0 


1.0 


1.35 


1.0 


Y 


dotted 


P/J 


DD 


(5) 


1.0 


1.0 


1.35 


1.5 


Y 


solid 



"NLO refers to the default next-to-leading order HVQMNR distributions. 
b LO refers to the default leading-order HVQMNR distributions. 
c cc refers to the default Pythia/Jetset cc distributions. 
d DD refers to the default Pythia/Jetset DD distributions. 
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Fig. 30. Sensitivity of single-charm and charm-pair distributions to higher- 
order perturbative corrections. The LO (dashed) and NLO (solid) distribu- 
tions are obtained from the HVQMNR generator; the parton-shower distribu- 
tions, with (solid) and without (dotted) intrinsic transverse momentum, from 
the Pythia/Jetset generator. See Table 5. 
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The next-to-leading order predictions and the parton shower prediction are 
also quite similar. The \Ap%\ and A<j> Pythia/ Jetset distributions with no in- 
trinsic transverse momentum included, indicate that the parton shower evolution 
is playing a very small role. The A<p parton shower distribution, in particular, is 
closer to the leading-order delta-function prediction than to the next-to-leading 
order prediction. Adding intrinsic transverse momentum, with Uk t = 0.44 GeV, 
brings the Pythia/Jetset prediction very close to the next-to- leading order 
HVQMNR prediction. 

A. 2 Sensitivity to the Mass of the Charm Quark 

In Fig. 31, we investigate the degree to which the single-charm and charm-pair 
distributions are sensitive to variations in the mass of the charm quark. All 
distributions are obtained from HVQMNR NLO calculations using the default 
values for all parameters — except for m c . Higher-order effects play a larger role 
as the charm-quark mass decreases because the ratio Q/Aqcd decreases, where 
Q is the energy scale of the interaction (Eq. 13). For the lightest charm-quark 
mass (m c = 1.2 GeV/c 2 ), the single-charm xp and p\ distributions are steepest 
because the outgoing charm quark can more easily radiate gluons; the single- 
charm y distribution is less peaked near y = 0; and the invariant mass of the 
charm-pair is significantly steeper than the higher mass predictions. The increase 
in higher-order effects for smaller m c is also evident in the A<p distribution, which 
is flattest for m c = 1.2 GcV/c 2 . 

A. 3 Sensitivity to Parton Distribution Functions 

In Fig. 32, we investigate the degree to which the single-charm and charm- 
pair distributions are sensitive to variations in the parton distribution functions 
and Aqcd- All distributions are obtained from HVQMNR NLO calculations 
using the default values for all parameters — except for the parton distribution 
functions. We examine predictions for four pairs of pion and nucleon parton 
distribution functions, sets (1) through (4) defined in Table 4. 

At fixed-target energies, the dominant contribution to the cc cross section is 
from gluon fusion. In Fig. 33, we compare the gluon distribution functions f g 
for sets (1) through (4). By energy conservation, the energy of the two colliding 
partons must be at least twice the mass of a charm quark to produce a cc pair; 
that is, > "^jf where \fs = 30.6 GeV is the center-of-mass energy of 

the colliding hadrons. Hence, for each set, the pion and nucleon functions are 

4 2 

obtained after imposing the constraint x^Xn > -^ JL . We impose this constraint 
because we want to investigate how the four sets compare in the region of x that 
we explore, not in the very low x region where the functions are largest. 

Although the four sets of parton distribution functions differ significantly, 
the single-charm and charm-pair distributions shown in Fig. 32 are not very 
sensitive to these differences. The sensitivity of the A<f> distribution is due to 
the variation in Aqcd in sets (1) through (4) (see Table 4). As the value of 
Aqcd increases, the ratio Q/ Aqcd decreases, where Q is the energy scale of the 
interaction (Eq. 13), causing higher-order effects to play a larger role. Hence, 
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Fig. 31. Sensitivity of single-charm and charm-pair distributions to variations 
in the mass of the charm quark. All distributions are obtained from HVQMNR 
NLO calculations with the default values for all parameters except m c . The 
values of m c used are 1.5 GeV (solid), 1.2 GeV (dashed) and 1.8 GeV (dotted). 
See Table 5. 
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Fig. 32. Sensitivity of single-charm and charm-pair distributions to variations in 
the parton distribution functions. Sets (1) (solid), (2) (dashed), (3) (dotted) and 
(4) (solid) are defined in Table 4. All distributions are obtained from HVQMNR 
NLO calculations with the default values for all parameters except the parton 
distribution functions. See Table 5. 
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PDF Set (1) 




Fig. 33. Comparison of the gluon distribution functions for sets (1) through (4), 
defined in Table 4. The pion and nucleon functions for each set are obtained from 
PDFLIB after imposing the constraint ^Jx^xn > "^jf , where ^/s is the center-of- 
mass energy of the colliding hadrons. The left plot shows the pion and nucleon 
gluon distribution functions for set (1). The middle and right plots show, for 
the pion and nucleon, respectively, the asymmetries between the set (1) gluon 
distribution function and each of the other three (n=2-4) gluon distribution 
functions. 

the flattest Acj> distribution results from using Set (4) {A { q ] cd = 300 GeV); the 
steepest A<p distribution, from using set (3) (Aq^ cd = 100 GeV). 

A. 4 Sensitivity to Factorization and Renormalization Scales 

In Fig. 34, we investigate the degree to which the single-charm and charm-pair 
distributions are sensitive to variations in the renormalization and factorization 
scales. All distributions are obtained using the HVQMNR NLO calculation. We 
set the two arbitrary scales equal to each other, /i = jip = hr, and obtain 
distributions for n = Q/2, Q, and 2Q, where Q gives the energy scale of the 
interaction (Eq. 13). We use the GRV parton distribution functions for both the 
pion and the nucleon (set (2) in Table 4), which have been evolved down to /i 2 , 
= 0.3 GeV 2 . With this choice, the factorization scale \i can go as low as m c /2 
without going below /z . As mentioned, the degree to which the distributions are 
sensitive to variations in the renormalization and factorization scales gives an 
indication of how much (or little) we can trust the o? a calculation. As expected, 
the distributions that are most sensitive to variations in \i are those distributions 
that are trivial at leading-order: |Z\p 2 |, A<j>, and p 2 -. The smaller the factor- 
ization and renormalization scales, the broader these distributions are. That is, 
the higher-order a 3 terms play a larger role, compared to the leading-order a 2 
terms, as renormalization and factorization scales decrease. The sensitivity of 
these scales indicates that the model may have several ways of obtaining ac- 
curate predictions, both by adjusting the model's parameters, and by adding 
higher order terms. 

A. 5 Sensitivity to Higher-Order Nonperturbative Effects 

In Fig. 35, we look separately at the effects of parton showers, intrinsic trans- 
verse momentum, and hadronization. All distributions are obtained using the 
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Fig. 34. Sensitivity of single-charm and charm-pair distributions to variations 
in the factorization and renormalization scales. All distributions are obtained 
from HVQMNR NLO calculations. Rather than using the default set of parton 
distribution functions, we use the GRV functions, which are evolved down to 
fj,Q = 0.3 GeV 2 . The values of fio/Q used are 0.5 (solid), 1.0 (dashed) and 1.5 
(dotted). The energy scale Q is defined in Equation 13. See Tables 4 and 5. 
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Pythia/Jetset event generator. The distributions obtained using the default 
settings (solid) include all three effects. We compare these default distributions 
to three sets of distributions that are obtained by including: 

• only hadronization, but no parton shower evolution or intrinsic transverse 
momentum (solid); 

• only the parton shower evolution, but no intrinsic transverse momentum 
or hadronization (dashed); 

• only intrinsic transverse momentum, but no hadronization or parton shower 
evolution (dotted). 

For the longitudinal momentum distributions (xf, Sxp 1 Axf, y, Sy, Ay), the 
most important factor is whether or not hadronization is included. The two 
sets of distributions that include hadronization effects are quite similar; the two 
sets of distributions that do not include hadronization effects are similar; but 
the latter two sets of distributions differ significantly from the former two sets. 
In the Pythia/Jetset hadronization model, the broadening of the longitudi- 
nal momentum distribution is the result of color-connecting the charm quark 
to a valence antiquark (or diquark) from one of the colliding hadrons and the 
anticharm quark to a valence quark from the other colliding hadron. 

All three higher-order effects broaden the leading-order delta function predic- 
tion for the A(j) distribution. The broadening due to the parton shower evolution, 
however, is significantly smaller than the broadening due to either the hadroniza- 
tion process or the addition of intrinsic transverse momentum (<Tfc t = 0.44 
GeV/c). The latter two effects broaden the Acj) distribution by roughly the same 
amount. 

All three higher-order effects also broaden the leading-order delta function 
prediction for the c - distribution. In this case, however, the broadening due 
to the parton shower evolution is larger than the broadening due to either 
hadronization effects or the addition of intrinsic transverse momentum {a^ t = 
0.44 GeV/c). 

A. 6 Sensitivity to Intrinsic Transverse Momentum 

In Fig. 36, we investigate the degree to which the single-charm and charm-pair 
distributions are sensitive to variations in the amount of intrinsic transverse 
momentum added to the hard-scattering partons that collide to form a cc pair. 
All distributions are obtained from the Pythia/Jetset event generator, with 
default settings for all parameters except for the width of the Gaussian intrinsic 
transverse momentum distribution, a kt ■ 

When intrinsic transverse momentum is included, the hard-scattering par- 
tons from the colliding hadrons arc no longer necessarily moving parallel to 
the colliding hadrons. The plane that is transverse to the axis of the parton- 
parton collision — which cannot be determined experimentally — is no longer 
the same as the plane that is transverse to the beam axis. Hence, including in- 
trinsic transverse momentum smears the leading-order prediction p t c = —p t c- 
Not surprisingly, the distributions that are most sensitive to variations in a kt 
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All Effects, Only PS, Only k,, Crly Hadr. 




Fig. 35. Sensitivity of the single-charm and charm-pair distributions to the 
parton shower evolution (PS), the addition of intrinsic transverse momen- 
tum, and the hadronization process. All distribution are obtained from the 
Pythia/Jetset event generator. The solid distributions include all three ef- 
fects; the dashed distributions include only the parton shower evolution; the 
dotted distributions include only intrinsic transverse momentum; the solid dis- 
tributions include only hadronization. See Table 5. 



http : //link. springer .de/link/ service/ journals/10105/ index.html 



EPJdiroct CI (2008) Springer- Verlag 



131 




Fig. 36. Sensitivity of single-charm and charm-pair distributions to variations in 
the amount of intrinsic transverse momentum added to the hard-scattering par- 
tons. All distributions are obtained from the Pythia/Jetset event generator, 
with default settings for all parameters except <Jk t - The values of <7k t used are 
0.44 GcV/c (solid), 0.7 GcV/c (dashed), 1.0 GcV/c (dotted) and 1.5 GcV/c 
(solid). See Table 5. 
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are those transverse distributions that are trivial at leading-order: A(f>, 
and vice- As the width <Jk t increases, these distributions become natter. 

A. 7 Summary 

In this section, we briefly summarize the results of the comparisons shown in 
Figs. 30-36. 

The longitudinal momentum distributions — xf, Sxf, Axf, y, Sy, and 
Ay — are relatively insensitive to all variations considered above, except for 
inclusion or omission of the Pythia/ Jetset hadroproduction hadronization 
(Fig. 35). The steepness of the invariant mass distribution is also sensitive to 
whether or not hadronization is included, as well as to the mass of the charm 
quark (Fig. 31). Therefore, the measured distributions for these physics variables, 
discussed in Sec. 6, provide a test of the Pythia/ Jetset hadronization model — 
in particular, a test of the string topology scheme that color-connects the charm 
quark to a valence quark from one of the colliding hadrons and the anticharm 
quark to a valence quark from the other colliding hadron. 

The transverse distributions Ep1, A(j>, and c - are sensitive to al- 

most all variations considered above because they are sensitive to the degree of 
correlation between the charm and anticharm transverse momenta. Varying m c 
(Fig. 31), Aqcd (Fig. 32), or hr (Fig. 34) in the next-to-leading order calcula- 
tion changes the definition of the running coupling constant a s , which is approx- 
imately proportional to 1/ \n(nn/ Aqcd)- As the coupling constant increases — 
that is, as m c decreases, Aqcd increases, or /x/j decreases — higher-order ef- 
fects play a larger role, and consequently the charm and anticharm transverse 
momenta become less correlated. The other methods we discussed for including 
higher-order effects were parton showers, intrinsic transverse momentum, and 
hadronization. 

In Sec. 6, we quantify the degree of correlation between the transverse mo- 
menta of the D and D mesons from our DD data sample. The sensitivity of 
the NLO predictions to the arbitrary renormalization and factorization scales 
(Fig. 34) indicates that higher-order perturbative corrections are important. In 
principle, one could determine the sets of theoretical parameters that generate 
predictions that are in good agreement with the full range of experimental re- 
sults. The set of fit values chosen, however, may not be unique. The fit values of 
(e.g., the mass of the charm quark) would depend on the values of the renormal- 
ization and factorization scales. For example, if a renormalization scale of Q/2, 
rather than Q, is assumed, then a smaller value for <Jk t or a larger value for m c 
could each be used to fit the data. 
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